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Chapter  I 


INTRODUCTION  | 

In  many  wave  propagation  problems,  diffraction  is  considered  to  be  a 
second-order  effect  that  diffuses  the  sharp  shadows  and  focal  surfaces  j 

of  ray  theory.  However,  if  the  frequency  is  low  enough  so  that  the 
wavelength  is  comparable  to  the  size  of  the  focusing  region,  structured 
fields  of  significant  levels  can  be  found  in  classical  shadows  and 
smooth,  featureless  enhancement  can  replace  sharp,  multiple  convergences 
common  at  higher  frequencies. 

This  investigation  examines  strong  diffraction  in  regions  of  complex 
refractive  focusing  by  quantitative  analysis  of  low-frequency  wave  pro¬ 
pagation.  Focal  surfaces  form  because  of  gradients  in  the  refractive 
index  such  as  are  common  in  the  atmosphere  or  the  ocean  (either  for 
acoustic  or  electromagnetic  waves)  and  the  frequency  range  considered  is 
such  that  the  wavelength  is  of  the  same  order  as  the  separation  of  these 
surfaces . 

Most  of  the  work  done  on  acoustic  (or  electromagnetic)  focusing  has 
been  in  simple  cases  where  normal-mode  expansion  (Ewing,  et  al.,  1957) 
is  practical  or  asymptotic  methods  (Sachs  and  Silbiger,  1971)  can  be 
used.  Long-wave  length  focusing  in  inhomogeneous  media  has  been  avoided 
for  good  reason:  if  the  wavelength  is  comparable  to  some  characteristic 
length,  the  wave  nature  of  the  field  dominates  and  neither  small  nor 
large  value  asymptotic  approximations  are  useful. 
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While  analysis  of  long-wavelength  focusing  can  be  difficult,  the  phe¬ 
nomenon  is  not  uncommon  either  in  ocean  acoustics  (Brekhovskhikh  and 
Lysanov,  1982)  or  atmospheric  acoustics  (Thomson,  1983).  Furthermore, 
when  conditions  exist  for  long-wavelength  focusing,  the  effects  of  per¬ 
turbations  on  signal  fluctuations  can  be  primarily  diffractive.  The 
usual  multiple  scattering  theories  assume  that  the  spatial  scale  of  the 
perturbations  is  either  much  smaller  or  much  larger  than  a  wavelength 
(for  example,  Skudrzyk,  1957  or  Flatte,  1979).  In  the  present  analysis, 
the  perturbation  scale  is  of  the  same  order  as  the  acoustic  wavelength 
and  so  these  limiting  cases  are  not  applicable.  Since  very  little  anal¬ 
ysis  of  strong  diffraction  effects  has  been  done  either  for  determinis¬ 
tic  focusing  or  for  perturbation -induced  fluctuations,  this  research 
considers  both  phenomena. 

Three  critical  areas  are  considered:  long-wavelength  focusing,  sig¬ 
nal  fluctuations  caused  by  perturbations  in  sound  speed,  and  the  tran¬ 
sition  in  fluctuation  strength  from  small  amplitude  variation  to  a 
large,  but  limiting,  amplitude  fluctuation.  Normally,  the  degree  of 
signal  fluctuation  is  related  to  the  sound-speed  variation  introduced  by 
the  perturbations.  There  comes  a  point,  however,  when  the  wavefront  has 
become  so  distorted  that  maximum  phase  interference  is  taking  place  and 
increased  perturbation  only  relocates  the  interference  minima  and  maxi¬ 
ma.  This  is  termed  the  saturation  regime  and  is  the  subject  of  the  last 
area  of  consideration. 

Deterministic  propagation  and  signal  fluctuation  are  usually  consid¬ 


ered  separately,  but  the  common  theme  of  strong  diffraction  suggests 
that  they  be  treated  together.  If  unperturbed  focusing  phenomena  were 
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the  only  interest  of  this  investigation,  the  problem  could  be  cast  in 
nondimens ional  terms.  Perturbations,  however,  are  peculiar  to  the  par¬ 
ticular  physical  problem.  To  treat  realistic  values  of  perturbation 
scale,  strength,  and  rate  of  change,  a  simple  physical  situation  is  mod¬ 
eled.  Thomson  (1983)  has  shown  that  strong  focusing  can  occur  in  the 
planetary  boundary  layer  and  Brown  (1980)  and  Thomson  have  provided  a 
detailed  set  of  meteorological  measurements  in  connection  with  noise  in¬ 
terference  studies  of  the  DOE/NASA  wind  turbine  generator  near  Boone, 
North  Carolina.  This  example  provides  an  excellent  basis  for  the  pres¬ 
ent  studies  in  long-wavelength  focusing  and  scattering  but  the  methodol¬ 
ogy  and  many  of  the  results  are  applicable  to  strong  diffraction  prob¬ 
lems  in  general.  Other  examples  include  long-range  refractive  focusing 
of  acoustic  waves  by  temperature  gradients  in  the  thermosphere,  short - 
and  long-range  focusing  of  sound  in  the  ocean,  and  ionospheric  refrac¬ 
tion  of  very  low  frequency  electromagnetic  radiation. 

Several  tools  were  developed  specifically  for  this  investigation. 
Much  of  the  analysis  of  the  unperturbed  field  is  performed  using  a  nor¬ 
mal-mode  solution  based  on  standard  theory  (for  example,  Ewing,  et  al., 
1957);  however,  because  as  many  as  4000  modes  must  be  calculated  in  some 
instances,  a  very  efficient  eigenvalue  search  procedure  was  required. 
Consequently,  a  matrix  estimation  scheme  introduced  by  Cooney,  et  al. 
(1981)  for  locating  bound-state  energy  levels  for  the  time-independent 
Schrodinger  equation  is  adapted  to  drive  an  iterative  search  procedure. 
The  resulting  technique  is  reliable  and  fast. 

Since  the  focal  surfaces  or  caustics  are  closely  spaced  in  terms  of 
the  acoustic  wavelength,  the  usual  first-order  asymptotic  correction  to 
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ray  theory  for  a  simple  caustic  produces  a  very  poor  approximation  to 
the  actual  field.  This  failure  prompted  the  development  of  a  second-or¬ 
der  correction  that  not  only  handles  the  multiple  focus  well  but  also 
provides  a  criterion  for  deciding  when  the  second-order  scheme  should  be 
used  and,  also,  when  the  detailed  structure  of  a  multiple  focus  will  be 
seen. 

In  order  to  solve  for  the  perturbed  field,  a  range -dependent  solution 
must  be  used.  This  requirement,  along  with  the  desire  to  investigate 
the  nature  of  the  field  over  a  broad  range  of  frequencies,  motivated  the 
development  of  a  finite-difference  approximation  for  the  time-dependent 
acoustic  wave  equation.  By  using  time-limited  pulses,  radiation  condi¬ 
tions  are  easily  satisfied  and  the  frequency  response  over  more  than  a 
decade  is  given  at  every  one  of  the  grid  points.  This  particular  tech¬ 
nique  introduces  no  artificial  dissipation  and  a  modification  suggested 
by  Vichnevetsky  and  Bowles  (1982)  makes  the  numerically-induced  disper¬ 
sion  controllable  and  virtually  independent  of  grid  orientation.  Of  the 
range-dependent  techniques  used  here,  the  finite-difference  method  is 
the  most  accurate  so  it  is  used  as  a  baseline  against  which  the  other 
methods  are  compared. 

Since  this  finite-difference  method  has  not  been  applied  previously 
to  time-dependent  acoustic  propagation  problems,  some  effort  is  spent  to 
verify  the  results.  During  this  validation  process,  a  new  solution  for 
time-dependent  wedge  diffraction  in  two  dimensions  was  developed.  The 
solution  is  based  on  the  method  of  normal  coordinates  outlined  by  Biot 
and  Tolstoy  (1957)  but  the  usual  sum  of  Legendre  functions  of  non-inte¬ 
ger  degree  has  such  poor  convergence  that  even  30  decimal-digit  preci- 
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sion  is  inadequate  to  calculate  it  except  in  very  special  cases.  Here, 
an  integral  transformation  is  introduced  and  the  result  is  a  simple  and 
rapidly  convergent  integral  of  elementary  functions  instead  of  a  summa¬ 
tion  of  special  functions.  This  wedge  solution  is  used  to  evaluate  the 
finite-difference  method  for  diffraction  into  a  shadow.  Several  other 
tests  are  performed  to  evaluate  the  boundary  conditions,  source  treat¬ 
ment,  and  refraction  treatment. 

Some  of  the  more  important  results  from  the  analysis  of  the  unper¬ 
turbed  field  concern  the  underlying  field  structure.  For  example,  there 
is  a  definite  interference  structure  in  the  shadow  zone.  Here,  diffrac¬ 
tive  spreading  from  classical  rays  on  the  illuminated  side  of  the  caus¬ 
tic  and  the  evanescent  "tails"  of  low-order  modes  combine  to  produce  a 
spatially  complicated  field.  In  contrast,  in  the  strong  enhancement  re¬ 
gion,  the  field  is  dominated  by  a  single  component  and  so  this  region  is 
generally  featureless.  The  transition  from  a  highly  structured  multiple 
focus  at  high  frequency  to  this  smooth  enhancement  region  at  low  fre¬ 
quency  is  illustrated  by  the  behavior  of  the  second-order  asymptotic 
correction  to  ray  theory.  Everywhere  in  the  field,  the  structure  is  de¬ 
fined  by  distinct  groups  of  modes  and  the  interference  of  these  groups. 

Once  perturbations  are  introduced  into  the  sound  speed,  the  signal 
levels  at  any  point  fluctuate  according  to  the  reaction  of  the  field 
components  to  these  perturbations.  While  the  fluctuations  resulting 
from  these  perturbations  can  be  calculated  easily  using  the  finite-dif¬ 
ference  method,  this  method  gives  no  physical  insight  into  the  specific 
origins  of  fluctuation;  therefore,  an  eigenvalue  shift  method  is  intro¬ 
duced.  Based  on  range-dependent  normal-mode  theory,  the  eigenvalue 
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shift  technique  isolates  the  major  effects  of  sound-speed  perturbation 
by  identifying  the  modal  components  of  the  field  that  are  most  affected 
by  the  perturbations.  Coupled  with  a  knowledge  of  the  underlying  field 
structure,  the  mechanisms  of  fluctuation  can  be  determined.  The  impor¬ 
tant  features  producing  the  fluctuations  are:  the  presence  of  a  dominant 
mode  group  identified  by  a  maximum  in  the  eigenvalue  shift  function;  the 
phase  shift  between  this  dominant  component  and  a  component  consisting 
of  the  remainder  of  the  mode  sum;  and,  the  relative  magnitudes  of  these 
two  components.  By  modeling  the  fluctuations  as  the  interference  of  two 
unequal  components  with  a  changing  relative  phase,  results  very  close  to 
those  of  the  complete  eigenvalue  shift  solution  are  obtained  with  sig¬ 
nificantly  less  computational  effort.  This  simplified  two-component 
model  gives  excellent  results  in  the  region  of  strong  enhancement  and 
weaker  but  adequate  agreement  with  the  finite-difference  results  in  the 
shadow  zone. 

Because  this  analysis  of  strong  diffraction  led  to  the  development  of 
a  number  of  tools,  the  actual  analysis  is  not  considered  until  Chapters 
V  through  VII.  The  first  several  chapters  develop  these  tools  and  es¬ 
tablish  their  validity.  If  only  a  general  overview  of  the  analysis  is 
desired,  Chapters  II  through  IV  can  be  skipped  on  the  first  reading  and 
referred  to  as  necessary  through  the  equation  references  in  the  subseq¬ 
uent  chapters. 


Chapter  II 


SOLUTIONS  FOR  THE  FIELD  NEAR  A  COMPLEX  FOCUS 

A  region  of  complex  focusing  in  which  there  are  several  areas  of  con¬ 
vergence  in  close  proximity  is  particularly  interesting  when  the  fre¬ 
quency  is  sufficiently  low  that  strong  diffraction  occurs.  In  this 
chapter,  the  basic  governing  equations  for  acoustics  are  briefly  re¬ 
viewed  along  with  a  simple  normal-mode  method  and  the  usual  first-order 
asymptotic  correction  to  ray  theory.  In  addition,  a  second-order  asymp¬ 
totic  correction  that  allows  calculation  in  the  vicinity  of  two  nearby 
caustics  is  developed.  The  tools  described  in  this  chapter  will  permit 
a  detailed  investigation  of  this  region  before  adding  the  complications 
of  perturbations. 

2 . 1  Governing  Equat ions 

The  equations  governing  the  behavior  of  acoustic  disturbances  in  a 
compressible  fluid  are  the  mass  conservation  equation 

tyj)  l  +  V-  (f,  vy  )  =  0  ,  (2.1) 

where  ^>r  is  the  density  and  v‘/  is  the  particle  velocity;  the  equation  of 
motion  for  an  inviscid  fluid  (neglecting  gravitational  effects) 
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<lvt  /} 


(V  V  )\ 


=  0 


(2.2) 


where  p,  is  the  total  pressure;  and  the  isentropic  relation  for  suffi¬ 
ciently  slow,  reversible  compression  and  rarefaction 


( ^  p/<^»  )s  =  c1 


(2.3) 


where  c  is  the  sound  speed. 

If  we  express  the  pressure,  density  and  particle  velocity  each  as  a 
static  term  plus  a  small  variation  (e.g.,  vj  =  vQ  +  v)  and  neglect  flow 
rotation  and  second-order  terms,  we  have 

^p/^t  +  ^  c*fv  +  v-^p  =  0  ,  (2.4) 


and 

+  7p/^#  +  V'Cv  v)  =  0  ,  (2.5) 

where  p  is  the  acoustic  pressure  and  v  is  the  acoustic  particle  veloci¬ 
ty. 

At  this  point,  we  can  neglect  the  bulk  flow  v0  (although  this  will  be 
reintroduced  in  Chapter  V),  differentiate  Equation  (2.4)  with  respect  to 


9 


time,  take  the  divergence  of  Equation  (2.5),  and  subtract  the  results  to 
get  the  usual  wave  equation  for  acoustics 


(2.6) 


It  is  interesting  to  note  the  consequences  of  this  last  manipulation. 
The  extra  differentiation  of  the  two  coupled  equations  in  p  and  v  is  a 
mathematical  artifice  designed  to  combine  them  into  a  single  equation. 
However,  this  new  equation  has  second  derivatives  in  both  space  and  time 
whereas  only  first  derivatives  appear  in  the  original  coupled  set.  Con¬ 
sequently,  discontinuities  in  first  derivatives  of  physical  properties 
are  not  admissible  solutions  to  this  new  equation  while  they  are  per¬ 
fectly  acceptable  in  Equations  (2.4)  and  (2.5). 

In  the  next  chapter,  we  will  briefly  consider  the  method  of  charac¬ 
teristics,  a  technique  applied  directly  to  Equations  (2.4)  and  (2.5) 
(or,  for  that  matter,  to  the  nonlinear  set  of  Equations  (2.1)  and 
(2.2)).  The  concept  of  propagating  discontinuities  in  the  first  deriva¬ 
tives  of  the  physical  properties  is  vital  to  the  development  of  signal 
flow  in  this  method.  Not  only  are  discontinuous  derivatives  allowed  but 
they  can  grow  to  be  so  large  that  the  properties  themselves  approach 
discontinuity  and  form  shock  fronts. 

Finally,  the  wave  equation  driven  by  a  point  source  at  za  is 


* 


£(>o  -  z0)  f(t) 


(2.7) 
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where  f(t)  is  the  source  driver  waveform. 


2 . 2  Norma l -mode  Methods 

2.2.1  Simple  two- layer  model 

The  mode  solution  in  two-dimensions  is  derived  by  taking  Fourier 
transforms  of  Equation  (2.7)  in  both  time  and  range  and  solving  the  re¬ 
sulting  inhomogeneous  differential  equation  by  variation  of  parameters. 
The  harmonic  pressure  amplitude  is  then  given  by 

Eu„(z>)  v„(z<) 

-  exp(iyx)  ,  (2.8) 

where  u  and  v  are  solutions  of  the  depth  equation 


r 

L 


d2/dz* 


=  0 


(2.9) 


and  u,  in  general,  satisfies  the  upper  boundary  condition  while  v  satis¬ 
fies  the  lower  condition.  The  Wronskian  WMir  of  u  and  v  is  zero  at  each 
eigenvalue  (At  an  eigenvalue,  u  and  v  each  satisfy  both  boundary 
conditions.)  The  source  and  receiver  depths  are  renamed  as  the  upper 
point  7.>  and  the  lower  point  z<.  In  the  cases  to  be  considered  here, 
all  of  the  modes  are  completely  trapped  and  there  is  no  branch  line  in¬ 
tegral  to  compute  (Brekhovskikh ,  1980).  Evanescent  modes  are,  however, 


important . 
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In  order  to  study  the  wave  structure  of  a  focused  field,  a  two-layer 
normal-mode  model  will  be  used  in  which  the  inverse  of  sound  speed 
squared  varies  linearly  with  depth  in  each  layer.  In  this  case,  the  so¬ 
lutions  to  the  depth  Equation  (2.9)  are  Airy  functions.  A  typical  two- 
layer  profile  is  shown  in  Figure  1.  The  physical  model  (i.e.,  the  model 
for  the  sound-speed  profile  of  the  atmospheric  boundary  layer)  intro¬ 
duced  in  Chapter  V  matches  this  profile  for  z/H  greater  than  zero  and 
less  than  unity.  Above  and  below  this  region,  the  sound  speed  is  con¬ 
stant  for  the  physical  model  (dashed  line)  while  the  profile  used  in  the 
normal-mode  model  (solid  line)  is  extended  until  the  sound  speed  becomes 
infinite.  The  gradient  of  the  lower  layer  is  very  small  so  the  errors 
introduced  by  this  mismatch  between  the  physical  model  and  the  normal¬ 
mode  model  are  negligible  except  for  ranges  well  beyond  those  considered 
in  this  study.  The  errors  introduced  by  the  mismatch  in  the  upper  layer 
are  confined  to  the  region  near  the  upper  boundary  (z/H  equal  to  unity) 
and  to  ranges  shorter  than  those  of  interest  here.  (For  the  purpose  of 
validating  the  finite-difference  method,  the  normal-mode  profile  is  mod¬ 
ified  in  Chapter  IV  by  terminating  the  layers  at  the  upper  and  lower 
boundaries  with  perfectly  reflecting  boundaries.) 

Since  the  sound  speed  goes  to  infinity  both  above  and  below  the  layer 
interface,  the  depth  functions  u  and  v  must  decay  exponentially  with 
distance  from  the  interface.  Consequently,  the  solutions  are 
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u(z)  =  Ai(Z,) 


v(z)  =  Ai(Z2) 


(2.10) 


where 


Z 


1,2 


( u,/c‘  - 


(2.11) 


and 


. 

(2.12) 

*/,l  =  (1/C*Z  '  1/C0  ,1  -  zo>  ■ 

The  eigenvalue  condition  is  given  by  setting  the  Wronskian  of  u  and  v  to 
zero,  or 


A  i  ( Z  # )  Ai(Zz)  -  ^Ai(Zz)  Ai(Z()  =  0  ,  (2.13) 


where  the  prime  denotes  the  derivative  with  respect  to  Z. 
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Once  the  eigenvalues  are  located  for  a  sufficient  number  of  modes, 
the  pressure  can  be  calculated  from  Equation  (2.8).  The  procedure  for 
computing  the  Airy  functions  is  discussed  in  detail  by  Gabrielson 
(1983). 

2.2.2  Matrix  iterat ion  method 

Because  of  the  severe  computational  load  imposed  by  using  many  sin¬ 
gle-frequency  mode  solutions  to  cover  the  spectrum  of  a  single  time-de¬ 
pendent  waveform,  a  very  efficient  mode-calculation  procedure  is  essen¬ 
tial.  Iterations  to  locate  exact  eigenvalues  are,  by  themselves, 
intensive,  requiring  many  calculations  of  Airy  functions  per  mode.  Fur¬ 
thermore,  finding  reasonably  accurate  estimates  to  initiate  the  itera¬ 
tion  process  can  be  just  as  time  consuming.  The  popular  WKB  method 
still  requires  complicated  calculations  and  can  be  effectively  used 
only  in  single  ducts. 

Cooney,  et  al.  (1981)  have  proposed  a  matrix  technique  for  estimating 
the  bound  states  of  the  one-dimensional  Schrodinger  equation.  This 
problem  is  identical  to  mode  location.  By  writing  simple  difference- 
equation  approximations  to  the  depth  Equation  (2.9)  at  many  uniformly 
spaced  depths,  a  tridiagonal  system  of  simultaneous  equations  is  formed. 
Many  efficient  computer  algorithms  exist  for  locating  the  eigenvalues  of 
such  systems  so  the  estimation  problem  is  rather  simple.  Even  if  the 
modes  are  degenerate  or  nearly  so,  the  matrix  technique  locates  all  ei¬ 
genvalues  in  order.  It  even  locates  evanescent -mode  eigenvalues  with 
equal  ease. 
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Once  the  estimates  have  been  computed,  it  is  a  simple  matter  in  a 
single-channel  profile  to  refine  these  estimates  by  Newtonian  iteration. 
The  Wronskian  derivative  is  used  to  correct  the  original  estimates;  this 
derivative  is  also  used  in  each  term  of  the  mode  sum  of  Equation  (2.8). 
The  entire  mode  solution  done  in  this  way  is  fast  enough  to  make  calcu¬ 
lation  of  thousands  of  modes  practical. 

Some  time  after  this  procedure  was  perfected,  Porter  and  Reiss  (1984) 
published  a  similar  method  for  ocean  acoustics.  They  estimate  the  ei¬ 
genvalues  in  the  same  way  but  iterate  to  refine  the  estimates  by  further 
matrix  operations.  If  the  sound-speed  profile  can  be  approximated  by 
some  tractable  functional  form,  Newtonian  iteration  is  considerably 
faster.  Even  for  arbitrary  profiles,  it  may  be  adviseable  to  use  a  cor¬ 
rection  scheme  like  that  of  Guthrie  (1974)  based  on  numerical  integra¬ 
tion  of  Equation  (2.9)  and  the  eigenfunction  orthogonality  condition. 

2 . 3  Asymptotic  Approx imat ions 

An  asymptotic  correction  to  ray  theory  is  described  by  Brekhovskikh 
(1980)  for  computing  the  field  in  the  vicinity  of  a  simple  caustic.  The 
development  of  this  correction  is  reviewed  below  in  order  to  provide  the 
framework  for  deriving  a  second-order  correction  for  computing  the  field 
in  the  vicinity  of  two  nearby  caustics.  We  will  follow  Brekhovskikh' s 
treatment  for  the  first-order  correction  except  that  we  will  use  an  ex¬ 
pression  for  the  field  in  two  dimensions  not  three.  This  only  involves 
dropping  a  factor  of  the  square  root  of  horizontal  distance. 

The  field  can  be  written  as  a  superposition  of  elementary  waves  as 


follows 
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so  that  at  the  saddle  point 


(2.16) 


Consequently,  one  way  of  finding  g  is  to  find  the  relation  between  x  and 
y  by  ray  tracing. 

A  representative  diagram  of  ray  range  as  a  function  of  horizontal 
wave  number  is  shown  in  Figure  2.  Where  the  slope  is  zero,  a  caustic 
forms:  there  are  two  in  this  diagram.  In  the  top  diagram,  the  region  at 
some  reference  wave  number  ^  is  approximated  by  a  straight  line  with 
the  appropriate  slope,  or 


Y*  Y 


Figure  2.  Typical  plot  of  range  as  a  function  of  wave  number  for  rays 
near  two  closely-spaced  caustics.  Dashed  lines  show  approximate 
curves  for  asymptotic  analysis:  from  top  to  bottom  they  represent 
simple  ray  theory,  first-order  caustic  correction,  and  second-order 
caustic  correction. 
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x  =  xo  +  x<»/(r  ■  Yo) 


(2-17) 


Here,  x'  is  the  derivative  of  x  with  respect  to  y  at  y.  Inserting 
this  equation  into  Equation  (2.16),  we  can  then  write  an  approximation 
to  the  field  as 


P  s:  exp(i^xfl) 


90 

/(^) /exP[i(x-x.)(r,)  -  ix'(^-y^)l/2j  dy  .  (2.18) 


If  we  deform  the  contour  by  rotating  it  by  -^'/4,  the  endpoints  will  co¬ 
incide  with  the  steepest  descent  path  and  the  integral  can  be  calculated 


—  (2 7r7±xo'  exp  £i(^x#jfl74  ±  Ql/4)  ,  (2.19) 


where 


Q  -  (x  -  x,)(2 /+x#' )’/2 


(2.20) 


The  sign  is  chosen  so  that  the  quantity  under  the  square  root  is  posi¬ 
tive  . 

It  is  instructive  at  this  point  to  define  a  ray  "width"  as  the  range 
difference  required  to  produce  a  TT/2  phase  change  in  the  field.  This 
can  be  evaluated  from  the  necessary  change  in  Q  in  Equation  (2.19) 
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(2.21) 


where  W  is  the  width  parameter.  Here,  W  gives  some  indication  of  the 
low-frequency  spreading  of  rays.  This  definition  is  limited  to  large 
values  of  the  derivative  should  really  be  evaluated  from  the  change 
in  coordinate  perpendicular  to  the  ray  path  rather  than  horizontal 
range. 

A  value  for  can  be  calculated  from  ray  theory.  For  x  equal  to 
xo,  the  field  intensity  is  proportional  to 


Pl  -  (2  ir/±x.J  )  £(jj) 


(2.22) 


The  field  can  also  be  calculated  by  standard  ray-tube  spreading  argu¬ 
ments 


.  (2.23) 


where  cs  and  c ^  are  the  sound  speeds  at  the  source  and  receiver  respec¬ 
tively  and  cv  is  the  vertex  speed.  Therefore 


KSir'k^c,  ^ CY/Cs  )*”  *  1  : 


.  (2.24) 


20 


This  expression  for  jjT  can  be  used  in  expressions  that  follow  to  relate 
the  asymptotic  expansions  for  the  field  integral  to  ray  theory. 

2.3.1  Simple  caustic 

At  a  caustic,  x^  goes  to  zero  and  the  field,  according  to  Equation 
(2.22),  becomes  infinite.  This  is,  of  course,  just  a  symptom  of  the 
breakdown  of  the  approximation  in  Equation  (2.17).  In  order  to  approxi¬ 
mate  the  field  near  a  caustic,  we  can  fit  the  range  curve  (Figure  2, 
middle  diagram)  with  a  parabola  at  the  caustic  point.  Then 

x  =  xa  +  x"  ^ I1  '  (2.25) 


21 


S 


Cx 


x„)(2/-x;) 


(2.28) 


Here,  the  real  root  of  the  cube  root  is  taken.  Now,  with  Equations 
(2.27)  and  (2.24),  we  can  compute  the  magnitude  of  the  field  near  the 
caustic. 

Let  us  consider  how  far  two  caustics  must  be  separated  for  them  to  be 
treated  as  distinct.  As  seen  in  Figure  2,  there  are  three  rays  at  any 
range  between  the  two  caustics  so  this  region  is  illuminated,  not  in 
shadow,  and,  from  the  shape  of  the  parabola  in  the  middle  diagram,  xe" 
is  negative  for  this  caustic.  The  zone  in  between  the  two  caustics  cov¬ 
ers  ranges  less  than  xc,  therefore,  from  Equation  (2.28),  S  is  negative 
there.  (S  is  also  negative  for  the  intervening  region  when  calculated 
for  the  other  caustic.)  For  negative  arguments,  Ai  is  oscillatory  and 
the  first  half-cycle  of  these  oscillations  is  complete  for  S  equals  -IK 
If  we  take  this,  somewhat  arbitrarily,  to  be  the  range  interval  CD  need¬ 
ed  for  the  caustic  to  be  distinct,  we  find  from  Equation  (2.28)  that 

’/j 

co  =  |  xa '/ 2  J 


(2.29) 
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2.3.2  Multiple  caustic 

If  there  are  two  caustics  closer  than  Cp,  then  a  further  modification 
to  this  asymptotic  theory  is  required.  Instead  of  fitting  one  caustic 
with  a  parabola,  we  can  fit  the  pair  with  a  cubic.  This  is  illustrated 
in  the  bottom  diagram  of  Figure  2.  The  range  equation  is 


x  =  x. 


+  x0'(r-fc)  +  xe"'(r-r#)3/6  . 


(2.30) 


and  the  field  integral  can  be  approximated  as 


(24/x"'  (^  )  exp(i^xo)  J exp[i(Yt-XtZ-tV)] 


dt  ,  (2.31) 


where 


(x."72  *)  (ffc)  . 


(2.32) 


and  the  parameters  X  and  Y  are  given  by 


*o'  (6/xJ") 


(2.33) 


(x  -  xQ)  (24/ x"' ) 
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The  integral  in  Equation  (2.31)  can  be  integrated  numerically  if  the 
contour  is  first  rotated  by  -ff/S.  The  convergence  is  then  quite  rapid. 

A  similar  integral  is  used  by  Pearcey  (1946)  to  evaluate  the  field 
near  a  cusp  which  is  just  a  limiting  case  of  two  caustics.  Pearcey  fits 
the  cubic  by  point,  slope,  and  first  and  second  derivative  at  the  cusp 
while,  in  the  above  analysis,  the  cubic  is  fit  by  point  and  slope  at 
each  caustic.  In  this  way,  the  caustic  locations  are  preserved;  how¬ 
ever,  they  are  forced  to  be  equal  in  strength. 

Incidentally,  a  plot  of  the  magnitude  of  the  integral  in  Equation 
(2.31)  is  given  by  Brekhovskikh  (1980).  The  integral  there  is  slightly 
different  (all  of  the  signs  in  the  exponent  are  positive),  but  the  re¬ 
sult  is  exactly  equal  to  the  complex  conjugate  of  that  obtained  here. 

These  asymptotic  approximations  will  be  compared  to  a  normal-mode  so¬ 
lution  in  Chapter  VI  where  we  will  see  that  the  multiple  caustic  approx¬ 
imation  does  generally  predict  the  correct  structure.  On  the  other 
hand,  the  first-order  asymptotic  correction  results  are  consistently 
poor. 

The  techniques  developed  in  this  chapter  will  be  applied  in  Chapter 
VI  to  the  investigation  of  the  field  near  a  complex  focus.  In  that 
chapter,  both  mode  theory  and  asymptotic  ray  theory  will  be  used  to  ana¬ 
lyze  the  structure  of  the  unperturbed  field.  The  next  chapter  (Chapter 
III)  will  summarize  the  tools  that  are  used  for  the  analysis  of  pertur¬ 


bations  . 


Chapter  III 


METHODS  FOR  RANGE -DEPENDENT  PERTURBATIONS 

3 . 1  Field  Solutions 

Once  spatial  perturbations  are  introduced  into  the  sound-speed  pro¬ 
file,  the  problem  becomes  range -dependent  and,  in  the  problem  considered  «■» 

in  Chapter  V,  strongly  so.  To  correctly  analyze  acoustic  propagation, 
the  solution  must  account  for  rather  large  variations  in  sound  speed 
over  one  wavelength.  We  will  first  consider  single- frequency  methods 
and  then  those  techniques  that  yield  a  direct  solution  in  the  time  do¬ 
main. 

3.1.1  Single-frequency  (CW)  techniques 

Ray  tracing  programs  have  been  developed  for  media  with  range  as  well 
as  depth  dependence  (Andersen  and  Kak,  1982);  however,  they  do  not  fully 
account  for  diffraction  even  with  wave  corrections  for  simple  caustics. 

Adiabatic-mode  formulations  (Pierce,  1965)  and  parabolic-equation  solu¬ 
tions  (Tappert,  1977)  can  cope  with  slow  variations  of  parameters  in 
range  but  they  are  s ingle- frequency  models;  broadband  representation  is 
difficult  and  expensive  (see  Chapter  IV).  In  addition,  phase  errors 
(for  example,  Brekhovskikh  and  Lysanov,  1982)  in  the  parabolic  method 
could  alter  the  physical  phase  fluctuations  near  the  focusing  regions. 

A  coupled-mode  approach  (Milder,  1969)  may  handle  the  variations  in 
range  encountered  but  only  at  an  even  greater  cost.  Note  that  range-de- 
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pendent  normal-mode  theory  will,  however,  be  used  to  estimate  fluctua¬ 
tions  . 

3.1.2  Time-dependent  solutions 

The  time-dependent  wave  equation  is  a  hyperbolic  partial  differential 
equation  and,  therefore,  can  be  reduced  by  the  method  of  characteristics 
to  a  simpler  problem  of  signal  progression  along  certain  surfaces  (von 
Mises,  1958).  In  one  space  dimension,  this  approach  works  very  well 
even  for  nonlinear  partial  differential  equations. 

As  an  example,  consider  a  plane  wave  incident  on  a  region  in  a  fluid 
in  which  the  sound  speed  changes  linearly  over  some  finite  distance  (see 
Figure  3).  There  is  no  impedance  discontinuity  to  produce  a  reflection 
but  there  is  weak  diffractive  backscatter  from  the  sound-speed  gradient. 
As  shown  in  the  figure,  the  method  of  characteristics  (dashed  line)  pro¬ 
duces  an  accurate  calculation  of  the  pressure  impulse  response  of  such  a 
layer  as  compared  to  a  much  more  expensive  superposition  of  single-fre¬ 
quency  plane  wave  reflection  responses  (solid  line).  Also  notice  that 
the  superposition  solution  suffers  from  oscillatory  "ringing"  at  the 
edges  of  sharp  level  changes  (Gibbs'  phenomenon). 

For  more  than  one  spatial  dimension,  Butler  (1960)  has  shown  that  the 
method  of  characteristics  can  still  be  used  but  the  resulting  algorithm 
introduces  a  great  deal  of  dissipation  (false  absorption)  and  dispersion 
(phase  speed  errors).  Furthermore,  the  code  is  complicated  and  boundary 
conditions  are  difficult  to  implement. 

Moretti  (1979)  has  developed  a  much  simpler  finite-difference  scheme 
that  is  weighted  so  that  signals  are  propagated  in  the  proper  character- 
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istic  directions;  the  method  has  proven  quite  successful  for  time-depen¬ 
dent  fluid  flow  problems.  For  acoustics,  however,  dissipation  in  the 
numerical  code  rapidly  distorts  the  waveshape.  This  will  be  illustrated 
below. 

One  important  advantage  of  the  method  of  characteristics  is  that  the 
physics  of  signal  propagation  is  respected;  the  numerical  grid  advances 
in  the  direction  of  local  signal  flow  rather  than  in  an  arbitrary  direc¬ 
tion.  Many  codes  advance  the  signal  using  data  from  field  points  that 
could  not  physically  affect  the  signal  (e.g.,  using  points  that  the  sig¬ 
nal  could  not  have  reached  yet).  This  turns  out  not  to  be  too  critical 
for  linear  acoustics  but  if  there  are  finite  amplitude  signals,  these 
will  be  seriously  distorted.  With  a  small  enough  grid  to  combat  the 
dissipation,  Moretti's  technique  may  be  useful  for  problems  in  nonlinear 
acoustics . 

Hyperbolic  partial  differential  equations  lend  themselves  to  analysis 
based  on  signal  dissipation  and  dispersion.  Numerical  solutions  can  be 
tested  by  introducing  the  Fourier  components  of  a  pulse  and  examining 
the  level  change  and  phase  shift  as  the  solution  proceeds  in  time  (von 
Neumann  and  Richtmyer,  1950).  Using  this  procedure,  a  number  of  finite- 
difference  solutions  to  either  the  wave  Equation  (2.6)  or  the  coupled 
governing  Equations  (2.4)  and  (2.5)  were  evaluated.  Among  them  were 
Euler  and  MacCormack  predictor-corrector  schemes  (MacCormack,  1970), 
first-  and  second-order  Moretti  algorithms  and  a  direct  finite-differ¬ 
ence  solution  of  the  wave  equation.  This  Fourier  analysis  is  outlined 
in  the  next  section  along  with  the  development  of  the  direct  solution 


method . 
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3.1.3  Direct  finite -dif ference  solutions 

The  selection  of  a  direct  finite-difference  approximation  to  the  wave 
equation  was  motivated  by  the  ability  of  this  approach  to  model  short 
duration  impulse-like  signals  (very  large  bandwidth)  in  a  range-depen¬ 
dent  environment  with  no  numerically  induced  dissipation  and  only  a 
small,  predictable  amount  of  dispersion.  In  addition,  for  low  frequen¬ 
cies,  the  method  has  about  the  same  computational  expense  as  a  normal¬ 
mode  solution  in  a  range- independent  duct.  This  latter  point  was  only 
evident  after  the  algorithm  was  fully  verified;  however,  the  dissipa¬ 
tion/dispersion  characteristics  can  be  determined  from  the  basic  approx¬ 
imation.  Hence,  we  will  consider  these  characteristics  along  with  the 
development  of  the  finite-difference  approximation. 

Some  special  notation  is  convenient  for  describing  finite-difference 
equations:  let  the  pressure  be  given  at  discrete  points  (in  space  and 
time)  by  the  symbol 

Pjm  =  P(l*x,  may,  nit)  .  (3.1) 


To  further  reduce  the  complexity,  denote  the  local  reference  point  as 
p°  =  p.  and,  for  example,  p'  =  p” 

The  wave  equation  in  two  dimensions  is 


W 


/  *• 

-  Ptt  / c 


+ 


(3.2) 
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By  replacing  the  partial  derivatives  with  their  first-order  finite-dif¬ 
ference  approximations  and  solving  for  the  pressure  at  the  (n+1)  time 
increment,  we  can  write 


2d  -  2qi)P;o  -  p^ 

+  +  Po,  +  Pi  +  O 


(3.3) 


where  q  =  c^t/^x  and  and  ay  are  equal  for  a  square  grid.  This  is  the 
fundamental  equation  for  computing  the  field  points  at  advancing  time 
but  there  are  several  ways  to  select  the  grid  points.  Figure  4  shows  a 
field  point  and  its  nearest  eight  neighbors.  If  the  x  axis  is  horizon¬ 
tal,  the  open  circles  are  used  as  grid  points  in  Equation  (3.3)  with 
ax  =  h.  If,  on  the  other  hand,  the  coordinate  system  is  rotated  45  de¬ 
grees,  the  dark  circles  are  used  with  ax  =  Jl  h.  Normally  dissipation/ 
dispersion  characteristics  are  functions  of  the  angle  of  propagation 
with  respect  to  the  grid;  however,  Vichnevetsky  and  Bowles  (1982)  have 
suggested  a  weighted  average  of  difference  equations  for  the  two  grids 
discussed  above.  The  resulting  equation  for  p  at  the  (n+1)  step  in  time 
is 


P~  = 

q  O  *  1  -  2* 

2Poo  -  Poo  +  q 

_ 

/2  + 

(i  y?)G0 

_  o  i  o  .  o 

=  p  +  p  +  p  + 

rO  1  r!0  r-/o 

O 

P„-,  - 

/  o 

4P.O 

(3.4) 

0,0  o 

—  p  +  p  +  D  + 

KJ-»  P-// 

0 

P-/-,  ' 

/  ° 

^Poo 

where  d  is  a  weighting  factor  that  ranges  from  0  to  1 . 
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Consider  the  effects  of  this  differential  equation  on  a  plane  wave 
propagating  through  the  grid.  Let 

Poo  =  exp[(o<  -ikc*)nat]  exp[i(kxl  +  k^rn)h]  ,  (3.5) 


where  k  is  the  wave  number  vector  with  components  kx  and  kj,  <*.  is  the 
apparent  absorption,  c*  is  the  apparent  phase  speed  and  h  is  the  grid 
spacing.  If  this  form  for  p  is  substituted  into  Equation  (3.4)  and 
solved  for  the  pressure  ratio  per  time  step,  the  wave  shape  change  is 
given  by 


P<Io /P^  =  exp  (©(a  t)  exp(-ikc#at)  .  (3.6) 


By  defining  an  amplification  factor 


(3.7) 


and  an  apparent  phase  speed 


c*  =  arg(p0'0  /p®0  )/k  a  t 


(3.8) 
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the  dissipation  and  dispersion  effects  can  be  separated.  Ideally,  the 
amplification  factor  A  should  be  unity  and  c**  should  be  equal  to  the 
true  local  sound  speed  c. 

For  Equation  (3.4) 


A 


D  + 


(3.9) 


and 


D 


(3.10) 


where  &  =  khcos  5/,  <5^  =  khsin^  and  </  is  the  angle  of  the  wave  number 
vector  with  respect  to  the  x  axis.  Here,  A  is  identically  unity  as  long 
as  the  magnitude  of  D  is  less  than  one.  As  described  below,  the  parame¬ 
ter  ^  is  selected  to  keep  the  dispersion  errors  independent  of  the  di¬ 
rection  of  propagation.  Also,  the  high-frequency  limit  of  the  source 
spectrum  sets  an  upper  limit  on  kh  for  a  particular  grid  spacing.  Con¬ 
sequently,  the  solution's  stability  is  controlled  by  q. 

For  a  hyperbolic  differential  equation,  the  Courant-Friedrichs-Lewy 
(Courant,  et  al.,  1956)  condition  (CFL  condition)  sets  an  upper  limit  of 
one  on  q.  (This  is  a  necessary  but  not  sufficient  condition  for  stabil- 
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ity.)  Physically,  the  CFL  condition  insures  that  the  numerical  solution 
keeps  up  with  the  signal  propagation.  The  parameter  q  is  the  ratio  of 
the  distance  the  actual  signal  travels  in  one  time  step  (cat)  to  the 
spatial  advance  of  the  grid  (h) .  If  the  grid  doesn't  advance  in  space 
faster  than  the  signal,  the  signal  will  eventually  pass  completely  out 
of  the  solution  region. 

In  practice,  q  must  be  somewhat  less  than  one  in  order  to  insure  that 
the  magnitude  of  D  in  Equation  (3.10)  remains  less  than  one.  A  practi¬ 
cal  limit  for  this  particular  finite-difference  scheme  is  0.8. 

The  apparent  phase  speed  ratio  is  given  by 


c*/c  =  tan-' (yj  1  -  D*'/D)/khq  (3.11) 

where  q  has  been  used  to  replace  the  &t  from  Equation  (3.8)  with  h. 
This  ratio  (which  should  be  unity)  is  plotted  in  Figure  5  as  a  function 
of  the  wavevector  angle  for  several  values  of  .  For  ^  of  roughly 
0.4,  the  dispersion  curve  is  very  nearly  isotropic.  In  Figure  6,  the 
value  =  0.4  is  taken  and  the  apparent  phase  speed  ratio  is  shown  as  a 
function  of  the  normalized  grid  size  kh/7r  for  various  values  of  q. 

This  direct  finite-difference  algorithm  introduces  no  artificial  ab¬ 
sorption  and  is  stable  as  long  as  the  time  step  is  chosen  to  restrict 
the  magnitude  of  q.  There  is  some  dispersion  but  it  can  be  limited  by 
proper  choice  of  the  spatial  grid  size  h  and  its  effects  are  virtually 
independent  of  the  direction  of  propagation. 


PHASE  SPEED  RATIO 

0.92  0.91  0.96  0.98  l.OC 


Figure  5.  Apparent  phase  speed  ratio  c*/c  as  a  function  of  the  angle  of 
the  wavefront  normal  with  respect  to  horizontal.  Each  curve 
corresponds  to  a  different  weight  ^  and  kh  =  ff/3. 
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By  way  of  comparison,  two  other  plots  are  shown  in  Figure  7  of  the 
dissipation  and  dispersion  characteristics  of  several  other  solution 
techniques.  A  simple  weighted-characteristics  code  is  shown  along  with 
a  simple  predictor-corrector  scheme  and,  for  reference,  the  curves  for 
the  direct  finite-difference  method.  Notice  that,  even  though  the  pre¬ 
dictor-corrector  scheme  is  a  second-order  method  (halving  the  step  size 
theoretically  cuts  the  error  by  a  factor  of  four),  its  dissipation  ren¬ 
ders  it  less  accurate  than  the  direct  method  which  is  only  first-order. 

3 . 2  Special  Methods  for  F luctuat ions 

Since  perturbations  of  an  acoustic  medium  such  as  the  atmosphere  or 
the  ocean  are  extremely  difficult  to  model  in  detail,  a  suitable  analyt¬ 
ical  method  for  determining  the  statistics  of  signal  fluctuations  based 
on  some  parameter  or  parameters  of  the  perturbations  would  be  most  use¬ 
ful.  Much  of  the  work  that  has  been  done  in  this  area  has  been  based  on 
models  of  homogeneous  and  isotropic  scatterers  (and  often  of  single 
scatterers)  in  an  isotropic  medium.  The  next  section  summarizes  some  of 
the  more  important  work  in  this  regard  but,  because  physical  perturba¬ 
tions  are  often  far  from  homogeneous  and  isotropic,  none  of  these  re¬ 
sults  are  directly  useful  here.  Section  3.2.2  summarizes  Flatte's 
(1979)  method  for  calculating  expected  fluctuation  statistics  for  an  in¬ 
homogeneous  medium  and  this  method  will  be  extended  to  regions  near 
caustics.  Finally,  a  somewhat  more  promising  approach  based  on  range- 
dependent  normal-mode  theory  is  developed. 
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3.2.1  Homogeneous  isotropic  turbulence  models 

Chernov  (1960)  and  Tatarskii  (1967)  formulated  the  fluctuation  theory 
for  single  spherical  scatterers  in  an  otherwise  homogeneous  and  unbound¬ 
ed  medium.  The  principal  result  was  the  identification  of  two  different 
dependencies  of  fluctuation  magnitude  on  range.  The  normalized  fluctua¬ 
tion  magnitude  V,  given  by 


V 


(3.12) 


'A 

varies  either  as  r  for  r  «  r  or  as  r  for  r  »  r  where  P  is  the 

O  & 

pressure  amplitude  of  a  particular  single-frequency  component 


r0  =  ka*/4 


(3.13) 


and  a  is  the  length  scale  of  an  individual  scatterer. 

By  analogy  to  Flatte's  (1979)  analysis  of  ocean  perturbations,  we  can 
identify  five  points  at  which  real  media  depart  substantially  from  the 
homogeneous  single-scattering  model: 

1.  multiple  scattering  may  predominate  in  real  turbulence; 

2.  the  medium  itself  introduces  refraction  thus  causing  determinis¬ 
tic  multipath; 

3.  the  strength  of  perturbations  varies  vertically; 

4.  because  of  turbulence  cascading,  the  dynamical  range  of  perturba¬ 
tions  is  large; 


5.  the  boundaries  introduce  multipath  by  reflection. 

Clearly,  the  single-scattering  model  cannot  predict  the  saturation  that 
is  observed  at  long  ranges:  V  approaches  a  limiting  value  of  0.52 
(Skudrzyk,  1957). 

Minor  improvements  have  been  made  to  this  theory:  for  example,  work 
on  absorption  by  scattering  by  Brown  and  Clifford  (1976)  and  the  work  on 
the  effects  of  turbulence  on  multipath  interference  by  Ingard  and  Maling 
(1963).  However,  these  models  correct  isolated  deficiencies  without 
treating  the  overall  complexity  of  wave  scattering  by  turbulence  in  a 
refracting  medium.  Roth  (1983)  demonstrates  the  inadequacy  of  these  and 
other  simple  models  when  applied  to  acoustic  fluctuations  in  an  atmos¬ 
pheric  boundary  layer  with  convective  instability. 

3.2.2  Strength  and  diffraction  parameters  for  f luctuat  ions 

One  method  for  analyzing  the  ability  of  a  realistic  medium  to  produce 
signal  fluctuations  has  been  developed  by  Flatte  (1979).  He  accounts 
for  sound-speed  perturbations  through  two  parameters:  a  fluctuation 
strength  parameter  &  and  a  fluctuation  diffraction  parameter  yl . 

The  strength  parameter  is  the  rms  phase  fluctuation  produced  by  per¬ 
turbations  along  the  unperturbed  ray  path.  The  diffraction  parameter,  a 
modification  of  Tatarskii's  (1971)  wave  parameter,  describes  how  much 
the  ray  path  must  be  changed  in  order  to  produce  phase  interference  at 
the  receiver  point.  Another  point  of  view  is  that  ^"contains  informa¬ 
tion  on  the  magnitude  of  the  perturbations  in  sound  speed  while  indi¬ 
cates  the  significance  of  multipath:  i.e.,  whether  distinct  paths  can 
form  or  the  single  path  is  dif fractively  "blurred". 
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This  formulation  accounts  for  the  inhomogeneity  of  the  medium  itself, 
for  the  distribution  of  perturbations,  for  multiple  scattering,  and  for 
multiple  paths.  Its  fundamental  defect  is  that  it  is  based  on  ray  theo¬ 
ry  and  cannot  easily  be  extended  into  classical  shadow  zones  or  to  low 
frequencies.  As  it  was  originally  developed,  it  cannot  be  used  near  or 
on  caustics,  but  this  deficiency  will  be  corrected  below. 

In  order  to  define  the  strength  parameter  $ ,  we  can  write  a  phase 
integral  over  the  ray  path  from  source  to  receiver  including  only  the 
phase  contributions  from  the  perturbation  component  of  the  sound  speed. 
The  mean  square  of  this  integral  averaged  over  time  is  the  fluctuation 
strength 

Z  r*‘i  Z 

§  =  <(k/^)>  .  O-l*) 

where  y*  is  the  normalized  perturbation  term 

/'‘(x.z.t)  =  £c(x,z,t)  -  cM(z)J  jca  ,  (3.15) 

and  cu  is  the  unperturbed  profile. 

The  diffraction  parameter  JY  is  considerably  more  difficult  to  com¬ 
pute  in  a  refracting  medium.  The  basic  definition  is 

(R/L)2/27T 


A 


(3.16) 
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where  L  is  the  vertical  correlation  distance  of  the  perturbation 
patches.  The  quantity  R  is  the  vertical  distance  through  which  an  in¬ 
termediate  point  along  the  unperturbed  ray  path  must  be  displaced  verti¬ 
cally  in  order  to  produce  a  TT  phase  difference  at  the  receiver.  In 
moving  this  intermediate  point,  the  source  angle  must  be  adjusted  so 
that  the  new  ray  hits  the  displaced  point  and  then  the  angle  must  be 
changed  at  that  point  so  that  the  ray  continues  on  to  the  receiver.  No¬ 
tice  that  the  strength  of  the  perturbations  does  not  affect  the  magni¬ 
tude  of  R;  it  is  purely  a  measure  of  diffraction.  The  local  yl  value 
is,  however,  weighted  by  the  normalized  local  $  and  integrated  along 
the  entire  path  to  get  the  composite  A.  .  This  prevents  a  very  large  ,/l 
at  a  point  at  which  there  are  very  weak  perturbations  from  biasing  the 
overall  /i  . 

In  practice,  the  path  deviation  required  to  produce  a  full  fr'  phase 
change  is  not  calculated  but,  instead,  the  intermediate  point  is  dis¬ 
placed  a  small  amount  and  R  is  computed  based  on  the  rate  of  phase 
change  with  vertical  displacement.  The  total  travel  time  along  a  ray 
can  be  expanded  as  a  series  for  small  intermediate  displacement  ^  as 
follows , 


/  x  //  3  ,,, 

T (£)  T(0)  +  £T  +  J  T  /2  +  £T  /6 


(3.17) 


The  eigenray  is  the  minimum  travel  time  ray  so  the  second  term  on  the 
right  side  is  zero  and  the  phase  difference  can  be  written  as 
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=  ^£2(t"  +  £T"73)/2 


(3.18) 


Normally  T  7  is  much  larger  than  T  for  small  £  but  T  vanishes  at 
caustics,  hence  the  third-order  term  will  be  retained.  When  the  phase 
difference  is  •Tf  then  is  equal  to  R;  therefore 


fR2(T"  +  RT"'/3)  =  1 


(3.19) 


gives  the  condition  for  R.  If  we  neglect  T  ,  the  expression  reduces 
to  Flatte's  expression  as  long  as  the  unperturbed  sound  speed  is  con¬ 
stant  everywhere.  (Flatte  derives  R  in  terms  of  path  length  and  that 
quantity  is  only  stationary  on  the  eigenray  if  the  eigenray  is  a 
straight  line.) 

At  this  point,  the  quantities  T<  and  T  must  be  related  to  the 
quantities  normally  obtained  by  ray  tracing.  Since  the  ray  vertex  speed 
c^  is  a  more  convenient  independent  variable  than  vertical  displacement, 
we  will  use  it  initially.  The  derivative  of  travel  time  T  with  respect 
to  c„  holding  x  constant  is  equal  to  the  derivative  of  T  holding  z  con¬ 
stant  (the  usual  ray  derivative)  minus  the  travel  time  increment  along 
the  extra  path  segment  of  the  latter  ray 


(<^T/^c^)x  =  ( <)  T/  9  cv)^  -  (dl/dx)(<3  x/  dcv)/cp 


(3.20) 
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where  cp  is  the  sound  speed  at  the  point  of  displacement.  By  Snell's 
law 


cp/cv  =  dx/dl 


dz/dx 


V 


(cv/cp)2 


1 


(3.21) 


so  that 


=  [< 


(  <?T/  )cv)  /(^.x/  <}c,) 


cv/c/ 


/V 


(Cy/cp)^ 


-  1 


(3.22) 


Now  T  is  identically  zero  on  the  proper  ray  between  the  source  and 
the  receiver  so  a  new  pair  of  rays  must  be  constructed  one  of  which  con¬ 
nects  the  source  to  a  point  offset  slightly  above  the  original  ray  and 
another  that  connects  that  offset  point  to  the  receiver.  This  is  re¬ 
peated  until  there  are  five  rays:  the  original  ray  and  four  other  rays 
that  have  vertical  midpoint  displacements  of  ‘p  ,  2^,  -^,  and  -2^. 
From  these  rays,  the  derivatives  T  "  and  T"  are  computed  numerically. 
Finally,  using  Equations  (3.19)  and  (3.16),  R  and  yl  are  calculated. 
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3.2.3  Range -dependent.  Normal-mode  Theory 

While  the  /[.  -  $  characterization  does  attempt  to  account  for  real 
distributions  of  sound-speed  perturbations,  it  is  based  on  ray  theory 
and,  as  such,  it  cannot  be  expected  to  produce  reliable  results  in  a 
highly  diffractive  field.  An  approximate  approach  can,  however,  be  de¬ 
veloped  from  a  wave  point  of  view  using  range-dependent  normal-mode 
theory . 

By  taking  the  Fourier  transform  in  time  of  the  wave  Equation  (2.7), 
the  Helmholtz  equation  for  a  medium  with  range-dependent  perturbations 
can  be  written 

d  p/  d  xA  +  +  kpZ  p  =  £(z-ze)£(x)  .  (3.23) 

Let 

p(x,z)  =  ^  %(*)  UJZ)  »  (3.24) 

m 

where  u^,  are  the  unperturbed  (i.e.,  range- independent)  eigenfunctions 
given  by  Equation  (2.9).  By  assuming  this  form  for  the  pressure,  mode 
coupling  is  neglected.  Substitute  Equation  (3.24)  into  Equation  (3.23), 
multiply  by  uw  and  integrate  over  z.  to  yield 

d2/^/dx2  +  (f^+  =  %(z,)  dz  *  (3.25) 


where 
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D 


2 


► 


dz 


(3.26) 


and  is  the  range-averaged  eigenvalue  used  to  compute  u„.  At  this 
point,  neglect  the  backscattered  field  by  writing  ^ as 


aw(x)  exp(iy^x) 


(3.27) 


If  the  problem  were  truly  range -independent,  the  would  be  constants. 
Substitute  this  form  for  ^  into  Equation  (3.25)  so  that,  after  neglect¬ 
ing  the  second  derivative  of  a„,  the  a^ functions  can  be  found  and 


W  exp 


dz 


(3.28) 


Thus,  the  principal  effect  of  the  perturbations  is  to  produce  a  shift  in 
the  eigenvalue  given  by 


Xr.f  '  Y„  *  /<<-  kl)%V2P"/“»dZ  '  <3':9> 

This  provides  at  least  an  approximate  method  of  computing  signal  fluctu¬ 
ations  based  on  wave  theory.  Both  the  eigenvalue  shift  technique  and 


the  extended  method  will  be  applied  in  Chapter  VII  to  a  specific  fo¬ 
cusing  problem. 


The  integrals  in  Equation  (3.29)  can,  of  course,  be  performed  numeri¬ 
cally  but  the  choice  of  Airy  functions  for  the  vertical  field  variations 
allows  these  integrals  to  be  evaluated  in  closed  form  for  either  con- 

«.  t- 

stant  or  linear  variation  in  k  -  kp  .  Bucker  (1970)  gives  an  exact  dif- 

L  L 

ferential  for  the  integrand  if  k  -  k p  is  constant 


uJdz 


-  d  (du/dZ)  - 


(3.30) 


If  we  use  this  expression  and  Stoke's  equation  for  the  Airy  functions 


d2u/dZ*  Zu  =  0 


(3.31) 


where  u  is  any  linear  combination  of  Ai  and  Bi,  we  can  write  an  exact 
differential  for 


3Zu2dZ 


(• 


d{u(du/dZ)  -  Z|^(du/dZ)  -  Zu 


(3.32) 


This  can  then  be  transformed  into  the  desired  form 


47 


i  l  , 
zu  dz 


u(du/dZ)  -  Z (du/dZ)  +  Z  u 


/3 


-  j^a/cf  -  y*)/yz  +  ?2o]  [(du/dz/-  zu2  j 


(3.33) 


Using  Equations  (3.30)  and  (3.32),  the  eigenvalue  shift  can  be  computed 
from  Equation  (3.29)  as  long  as  the  perturbations  can  be  approximated  by 

t  i 


segments  of  constant  and  linear  k  -k^  .  For  small  perturbations,  the 
relationship 


»  2<o*£c/c3 


(3.34) 


holds . 

In  the  next  chapter,  extensive  tests  of  the  finite-difference  method 
will  be  outlined  in  order  to  establish  the  technique's  validity.  The 
techniques  described  above  (the  normal-mode  and  the  finite-difference 
techniques,  in  particular)  will  be  applied  in  Chapter  VII  to  the  study 
of  fluctuations  induced  by  turbulence  in  a  region  of  strong  focusing. 


Chapter  IV 


VALIDATION  OF  FINITE-DIFFERENCE  SOLUTION 

Since  direct  finite-difference  methods  have  not  generally  been  ap¬ 
plied  to  time-dependent  acoustic  fields,  the  technique  should  be  tested 
for  its  accuracy  in  relation  to  other  known  solutions.  While  the  dis¬ 
persion/dissipation  analysis  of  the  previous  chapter  sets  some  specific 
limits  on  the  distortion  introduced  by  the  finite-difference  approxima¬ 
tion,  comparisons  against  known  solutions  are  necessary  to  insure  that 
diffraction,  refraction,  boundary  effects,  and  the  source  term  are  prop¬ 
erly  represented. 

4 . 1  Boundary  Conditions  and  Source  Diffusion 

Transparent-boundary  conditions  (Appendix  A)  are  used  in  the  finite- 
difference  model  to  avoid  contamination  of  the  focused  field  by  simply 
reflected  arrivals.  Consequently,  a  very  simple  test  can  be  performed 
by  keeping  the  sound  speed  constant  and  by  applying  transparent-boundary 
conditions  at  the  horizontal  boundaries  indicated  in  Figure  1.  The  fi¬ 
nite-difference  solution  can  then  be  compared  to  two-dimensional  free- 
field  spreading.  This  will  test  the  model's  representation  of  the 
transparent -boundary  conditions  and  of  the  diffusive  "tails"  associated 
with  two-dimensional  geometries. 

The  governing  equation  throughout  this  study  is  the  forced  wave  Equa¬ 
tion  (2.7).  In  this  first  test,  the  sound  speed  is  constant  everywhere 
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and  there  are  no  boundaries  so  by  converting  the  equation  to  cylindrical 
coordinates  with  the  source  at  the  origin  and  taking  the  Fourier  trans¬ 
form  in  time  and  the  Hankel  transform  in  space,  we  can  write 


P(4t»)  =  iF(^j)  (kx)  /  4 


(4.1) 


where  P  is  the  Fourier  transform  of  pressure  and  F  is  the  Fourier  trans¬ 
form  of  the  source  waveform 


f(t)  =  A  exp(-b2 12/4) 


(4.2) 


or 


F(*0 


(4.3) 


The  factors  A  and  b  scale  the  amplitude  and  the  extent  of  the  spectrum 
respectively. 

The  relative  gain  (or  loss)  of  the  finite-difference  model  with  re¬ 
spect  to  the  free-field  spreading  Equation  (4.1)  is  shown  in  Figure  8 
for  several  observation  points  roughly  midway  between  the  two  boundaries 
and  at  ranges  of  2.5  to  3.5  times  the  boundary  separation.  The  indepen¬ 
dent  variable  here  is  the  normalized  grid  spacing  kh.  The  steep  depar¬ 
ture  below  kh  equal  to  0.075  is  caused  by  the  breakdown  of  the  transpar- 


Figure  8.  Level  predicted  by  finite-difference  calculation  relative  to 
free-field  spreading  for  a  constant  sound-speed  channel  with 
transparent  boundaries.  The  *  marks  the  point  at  which  the  channel 
width  is  half  the  wavelength. 
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* 


ent-boimdary  condition:  below  this  point  the  boundaries  are  less  than  a 
half-wavelength  apart  and  errors  in  the  boundary  conditions  become  large 
enough  to  affect  the  solution  significantly.  Between  0.075  and  0.95 
there  is  a  small  departure  from  the  desired  response.  These  changes,  as 
long  as  they  are  smooth,  are  of  no  concern  for  this  study  since  enhance¬ 
ment  near  the  caustic  and  loss  in  the  shadow  will  cause  far  greater 
swings  in  magnitude  than  these  errors. 

4 . 2  Corner  Diffraction 

Although  the  physical  model  discussed  in  Chapter  V  does  not  have  any 
wedge-shaped  boundaries,  the  ability  of  the  technique  to  predict  the 
field  in  the  shadow  behind  a  wedge  is  an  important  indication  of  its 
overall  representation  of  diffraction.  This  is,  in  fact,  a  good  way  of 
separating  diffraction  and  refraction  since  this  shadowing  occurs  even 
for  constant  sound  speed. 

4.2.1  Integral  solution  for  diffraction  by  a  wedge 

An  analytical  solution  for  time-dependent  two-dimensional  wedge  dif¬ 
fraction  can  be  constructed  by  the  method  of  normal  coordinates  (Biot 
and  Tolstoy,  1957).  For  the  polar  coordinate  system  shown  in  Figure  9 
with  the  origin  at  the  apex  of  a  pressure-release  wedge,  the  solution  to 
the  forced  wave  equation  is 
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P(r,@,t)  =  c/^wV/r  r0'  • 

OO 

2P^-f(A)  Sin7/A  sinz/^ 

Vm  / 


(4.4) 


for  ct  greater  than  (r  -  ra)  and  less  than  (r  +  re)  and 


p(r,  G  ,t)  =  -2c lv9wy/tro  • 


OO 


sin  v„9 


(4.5) 


for  ct  greater  than  (r  +  ra)  where  a  unit  impulse  driver  has  been  used 

(i.e.,  F  =  1).  The  pressure  is  zero  before  the  range  of  Equation  (4.4). 
Also 


A 


(r* 


ro 


c2t2)/2rr„ 


(4.6) 


and 


n  rf  1 9 


(4.7) 
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This  is  the  exact  solution  for  a  pressure-release  wedge.  However, 
computations  of  p  are  not  practical  because  of  the  very  poor  convergence 
of  the  series  coupled  with  the  difficulty  of  computing  the  Legendre 
functions  for  non-integer  degree. 

This  problem  can  be  circumvented  by  transforming  the  expression  as 
follows.  First,  the  Legendre  functions  are  replaced  by  contour  inte¬ 
grals  of  simpler  functions  and  then  the  summation  and  the  integration 
are  interchanged.  This  allows  closed  form  evaluation  of  the  sum.  A 
simple  integral  results  that  is  easily  evaluated  numerically. 

Consider  Equation  (4.4):  Gradshteyn  and  Ryzhik  (1965)  give  an  inte¬ 

gral  form  for  the  Legendre  function  of  the  first  kind 

Tf 

Pv»*r(A)  =  -'/A  +  Va1-  l’  cos?)*'*  d£  ,  (4.8) 

o 

which  is  valid  for  Re(A)  greater  than  zero.  The  limit  on  ct  in  Equation 
(4.4)  corresponds  to  A  between  1  and  -1  so  Equation  (4.8)  is  not  com¬ 
pletely  suitable. 

Equation  (4.8)  is  limited  because,  for  Re(A)  less  than  zero,  the 
branch  cut  associated  with  the  branch  given  by 

A  +  yjA1-  -  l'  cos/  =  0  ,  (4.9) 

intersects  the  original  integration  path.  By  defining  a  new  integration 
contour  Cp  that  goes  from  0  to  tt/2,  then  from  Tf/l  to  w/ 2  +  i^  along  one 
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side  of  the  branch,  around  the  branch  point  (a  zero  not  a  pole),  down 
the  other  side,  and  finally  along  the  real  axis  to  77?  we  can  write 


(4. 10) 


for  the  entire  range  required  by  Equation  (4.4).  The  imaginary  part  of 
£  at  the  branch  point  is  given  by 


;inh 


(4.11) 


For  Re(A)  greater  than  zero,  this  diversion  is,  of  course,  unnecessary. 

Next,  substitute  this  integral  form  for  the  Legendre  function  into 
Equation  (4.4)  and  interchange  the  integration  and  summation  to  write 


V 


(A  + 


i 


cos 


(4.13) 


If  we  expand  the  sine  functions  as  exponentials  and  replace  V"  by 
exp(n*lnV),  then  the  entire  quantity  within  the  summation  can  be  written 
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as  terms  of  the  form  exp(n*())  and  the  sum  can  be  evaluated.  The  re¬ 
sulting  expression  for  the  pressure  is 

P  =  C  /V'  ^  (T_  -  T+)  d£ .  (4.14) 

C, 


where 


T+ 


(cos/i 


2  V  cos^ 


+  V  ) 


> 


(4.15) 


=  rr{0  +  9e)/&w 


To  transform  Equation  (4.5),  an  integral  representation  for  the 
Legendre  function  of  the  second  kind  is  used 

©O 

("A)  =  yi-  A  +  yjk2  -  l’cosh^)*  S  d£  .  (4.16) 

o 

For  the  region  of  validity  of  Equation  (4.5),  there  are  no  problems  with 
branch  cuts  except  for  selecting  the  overall  sign  of  the  -v^-1/2  root. 
The  easiest  way  to  pick  the  proper  sign  is  by  analytic  continuation, 
which,  in  this  case,  amounts  to  matching  the  two  pressure  expressions  at 
ct  =  r  +  re .  The  negative  root  is  the  proper  one. 
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Otherwise,  the  transformation  is  identical  to  the  previous  one  and 
the  result  is 


where 


c!2Tr&„s/rc^  • 


OO 

/(V  )  (T  _  +  T_+  -  T+_  -  T^+)  d£ 


(4.17) 


V 


/ 


(-  A  +  VAZ  -  1  cosh^  ) 


(4.18) 


Here,  T  is  identical  in  form  to  Equation  (4.15)  with 


Xj  =  rr  {0±Oo±rr)!9w 


The  solution  is  now  formally  complete;  however,  the  computation  of 
the  infinite  integral  in  Equation  (4.17)  is  easier  in  two  stages. 
First,  the  integration  is  performed  numerically  from  0  to  some  large  ^ . 
Then,  the  remaining  section  is  evaluated  by  approximating  the  integrand 
for  large 
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r 


(4.20) 


This  integration  can  be  performed  with  the  result  that 


2c  si  n(Tr&/&„)  sin  (n9J9w)  exp  [~ &(<?„+ 1;)]  cos  {-rrz/0„) 

«ra.<£.j)v^r  [(A*-i)/4]’r/J'’-"/4' 


(4.21) 


A  practical  choice  for  is  10  and  the  step  size  for  numerical  integra¬ 
tion  from  0  to  can  be  surprisingly  large  (0.5  was  used  with  simple 
trapezoidal  integration  in  the  example  below). 

4.2.2  Comparison  with  finite-difference  mode  1 

In  Figure  10,  the  time-domain  waveform  calculated  by  the  finite-dif¬ 
ference  method  is  compared  to  that  calculated  by  the  integral  solution 
for  a  3ff72  pressure-release  wedge.  The  receiver  point  is  well  within 
the  shadow.  Since  Equations  (4.14)  and  (4.17)  assume  an  impulse  source, 
the  results  from  these  calculations  have  been  convolved  with  the  source 
driver  waveform.  Equation  (4.2),  used  in  the  finite-difference  method. 


0.00  0.02  0.04  0.06  0.08  0.10 
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Ct 


Figure  10.  Diffraction  into  the  shadow  of  a  pressure-release  wedge. 
Finite-difference  calculation  is  shown  as  a  solid  line  while  the 
boxes  give  several  points  from  the  integral  solution.  Inset  gives 
the  x  and  y  coordinates  of  the  source  and  receiver. 
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Several  other  receiver  locations  were  evaluated,  all  with  equally 
good  results.  As  noted  before,  both  long  wavelength  diffraction  into 
shadows  and  two-dimensional  diffusion  must  be  handled  accurately  for 
these  results  to  agree  so  well. 

4.3  Focusing  at  Simple  and  Cusped  Caustics 

Up  to  this  point,  the  validation  tests  have  centered  on  those  wave 
effects  not  related  directly  to  refraction.  Since  the  region  of  inter¬ 
est  in  this  study  is  one  of  strong  refractive  focusing,  a  test  of  the 
finite-difference  method  under  these  conditions  is  essential.  An  appro¬ 
priate  comparison  solution  for  this  task  is  a  normal-mode  solution. 
Some  of  the  problems  in  using  a  superposition  of  single-frequency  solu¬ 
tions  to  check  a  time-domain  solution  will  become  evident  in  this  sec¬ 
tion.  Nevertheless,  reasonably  successful  comparisons  are  possible. 

For  these  comparisons,  the  sound-speed  profile  shown  in  Figure  1  was 
used  and  perfect,  pressure-release  boundary  conditions  were  applied  at 
the  upper  and  lower  boundaries  in  both  the  finite-difference  model  and 
the  normal-mode  model.  Since  the  boundaries  are  perfectly  reflecting, 
the  normal-mode  solution  predicts  infinite-Q  resonances;  these  would 
not,  of  course,  be  described  by  any  time-limited  solution.  In  order  to 
avoid  this  mismatch,  a  very  small  amount  of  attenuation  was  added  to  the 
modes  with  ray-equivalent  angles  greater  than  85  degrees  from  the  hori¬ 
zontal.  Also,  because  the  driver  is  a  transient  and  the  range  is  rela¬ 
tively  short,  low-frequency  evanescent  modes  were  included  in  the  solu¬ 


tion. 
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Figure  11  shows  the  predictions  of  the  mode  model  (dashed  lines)  and 
the  finite-difference  model  (solid  lines)  at  a  cusp  (upper  plot)  and  on 
a  simple  caustic  (lower  plot).  Overall,  the  agreement  shown  is  excel¬ 
lent.  Notice  that,  at  the  upper  end  of  the  spectrum,  the  curves  are 
shifted  slightly.  This  is  evidence  of  the  phase  errors  introduced  by 
numerical  dispersion  in  the  finite-difference  method.  More  than  4000 
modes  are  included  in  the  normal -mode  solutions  shown  in  these  figures. 
If  the  evanescent  modes  are  neglected,  the  spectra  differ  markedly. 
Also,  without  the  attenuation  applied  to  the  higher  order  modes,  very 
large  peaks  appear  at  several  frequencies  in  the  mode  curve. 

The  efforts  of  this  section  underscore  some  of  the  problems  associat¬ 
ed  with  computing  time-dependent  solutions  by  Fourier  superposition  of 
single-frequency  solutions.  Two  fundamental  problems  limit  the  accuracy 
of  such  calculation.  The  first  is  the  reciprocal  relation  of  frequency 
resolution  and  signal  time  duration.  Many,  closely  spaced  single-fre¬ 
quency  solutions  are  required  if  a  long  time  duration  is  expected  in  the 
signal . 

The  second  problem  has  its  roots  in  the  type  of  equation  being 
solved.  The  single-frequency  solution  is  a  solution  to  an  elliptic  dif¬ 
ferential  equation  (the  Helmholtz  equation)  and,  therefore,  the  entire 
space  must  be  considered  out  to  all  boundaries  including  those  at  infin¬ 
ity.  A  change  anywhere  in  the  space  changes  the  solution  so  some  care 
must  be  exercised  in  selecting  boundary  conditions  even  if  the  real 
time-dependent  signal  would  not  reach  those  boundaries  during  the  obser¬ 
vation  period.  The  time-dependent  solution  is  the  solution  of  a  hyper¬ 
bolic  equation  (the  wave  equation)  and  its  progression  follows  the  ad- 


Lgure  11.  Comparison  of  finite-difference  (solid  lines)  and  normal- 
mode  (dashed  lines)  solutions.  Upper  plot  shows  the  field  level 
(arbitrary  reference)  at  a  cusp  while  lower  plot  gives  the  level  on 
a  simple  caustic. 
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vance  of  the  physical  wavefront:  conditions  outside  this  envelope  do  not 
affect  the  solution  in  any  way.  This  conflict  shows  up  in  the  apparent¬ 
ly  perfect  resonances  in  the  mode  solution  for  ideal  boundaries. 

* 

4.4  Limitations  of  the  Finite-difference  Solution 

k  -----  -  -  -  . " 

|  The  most  serious  limitation  of  the  finite-difference  technique  is  the 

dispersion  that  forces  a  relationship  between  step  size  and  maximum  fre¬ 
quency.  At  about  five  points  per  wavelength  (about  30Hz  in  the  atmos- 

t  pheric  problem  considered  later),  the  distortion  becomes  noticeable. 

I 

This  should  be  true  in  any  application  of  the  method  so,  given  the  maxi¬ 
mum  frequency  and  the  minimum  sound  speed  (to  get  the  minimum  wave¬ 
length)  the  grid  size  is  fixed.  In  pushing  the  limiting  frequency  up¬ 
wards,  be  aware  that  the  running  time  increases  as  the  linear  density  of 
grid  points  cubed  since  the  time  step  must  also  be  decreased  as  the  grid 
spacing  decreases. 

The  low-frequency  limit  is  forced  by  imperfections  in  the  transpar¬ 
ent-boundary  condition.  As  shown  in  Figure  8,  this  limit  is  kh  roughly 
equal  to  0.075  (about  2Hz  for  the  atmosphere  problem).  If  another  type 
of  boundary  such  as  a  rigid  or  pressure-release  boundary  is  acceptable, 
this  limit  would  be  much  lower.  It  could  also  be  lowered  by  moving  the 
boundary  farther  away  at  some  increase  in  cost:  doubling  the  distance 

to  the  boundary  would  only  double  the  computation  time. 

Given  these  limitations  (primarily  economic),  the  finite-difference 
scheme  is  well  suited  for  impulse  studies  in  two-dimensions.  Wave  ef¬ 
fects  are  accurately  represented  and  large  range  variations  as  well  as 
depth  variations  in  sound  speed  and  boundary  geometry  are  easily  accomo¬ 
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Chapter  V 


SOLUTION  CONTEXT 

Focusing  of  waves  occurs  frequently  in  nature  and,  in  a  number  of 
these  cases,  the  relevant  frequencies  can  be  in  the  range  of  strong  dif¬ 
fraction.  For  acoustic  waves  in  the  atmosphere,  focusing  can  take  place 
over  distances  ranging  from  less  than  a  hundred  meters  to  thousands  of 
kilometers.  Intense  temperature  inversions  over  snow  coupled  with  wind 
shear  can  produce  dramatic  short  range  focusing  (Thomson,  1984).  In 
this  case,  the  caustics  would  only  be  separated  by  a  few  meters  and 
strong  diffraction  could  extend  to  several  hundred  hertz.  At  he  other 
extreme,  acoustic  waves  from  the  SST  (Weber  and  Donn,  1982)  or  meteors 
(Bedard  and  Greene,  1981)  can  be  focused  by  strong  temperature  gradients 
in  the  thermosphere  although  here  the  focusing  region  may  be  on  the  or¬ 
der  of  1000km  from  the  soiree  and  strong  diffraction  would  only  be  ob¬ 
served  for  extremely  low  frequencies  (well  below  1Hz).  Of  intermediate 
scale  is  focusing  by  wind  shear  in  the  atmospheric  or  planetary  boundary 
layer  (PBL) . 

Also,  short  and  long  range  focusing  in  the  ocean  results  from,  for 
example,  surface  ducts  or  the  SOFAR  channel  (Urick,  1979).  Oceanic 
sound  speed  is  four  to  five  times  higher  than  sound  speed  in  the  atmos¬ 
phere  so,  for  the  same  range  scales,  the  upper  limit  in  frequency  for 
strong  diffraction  in  the  ocean  is  four  to  five  times  higher  in  the 
ocean  than  in  the  atmosphere. 
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At  very  low  frequencies,  ionospheric  refraction  of  electromagnetic 
radiation  (Brekhovskikh,  1980)  can  be  strongly  diffractive.  At  very 
high  frequencies,  refraction  and  focusing  that  create  optical  mirages 
(Fraser  and  Mach,  1976)  would  normally  be  well  represented  by  simple  as¬ 
ymptotic  correction  to  ray  theory. 

5 . 1  Focusing  in  the  Planetary  Boundary  Layer 

While  focusing  in  a  stratified  medium  could  be  analyzed  in  a  general 
way  with  normalized  distances  and  frequencies,  the  introduction  of  real¬ 
istic  levels  of  perturbations  depends  on  the  specific  environment  being 
studied.  In  order  to  consider  a  useful  perturbation  problem  and  to  be 
able  to  assign  values  to  parameters,  a  specific  problem  is  examined 
here:  short-range  focusing  in  the  PBL.  All  of  the  general  results  are, 
however,  still  applicable  to  any  strongly  diffractive  focusing  environ¬ 
ment  . 

5.1.1  Effects  of  wind  on  sound  speed 

Because  bulk  flow  (wind)  in  the  atmosphere  is  not  generally  negligi¬ 
ble  compared  to  the  speed  of  sound,  the  wave  equation  must  account  for 
these  flow  effects.  We  could  combine  Equations  (2.4)  and  (2.5)  into  a 
wave  equation  without  dropping  the  bulk  flow  terms  but  it  will  be  easier 
to  include  a  simplifying  assumption  if  we  instead  compute  the  implied 
dispersion  relationship.  By  substituting  time-harmonic  expressions  for 
p  and  v  (through  a  velocity  potential),  we  can  find  the  condition  that 
allows  non-trivial  solution  of  the  coupled  Equations  (2.4)  and  (2.5) 
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u)  IV.  =  c  +  v0-k/k  ,  (5-1) 

where  k  is  the  wave  number  vector  and  oo  is  the  angular  frequency.  In 
other  words.  Equation  (5.1)  gives  an  effective  sound  speed. 

In  the  atmosphere,  the  bulk  flow  is  primarily  horizontal  and  the  rays 
are  nearly  horizontal  particularly  at  the  altitudes  of  the  higher  wind 
speeds.  Thus,  we  can  write  an  approximate  relation  for  the  effective 
sound  speed 

=  c  +  U  ,  (5-2) 

where  U  is  the  component  of  the  wind  in  the  direction  of  propagation. 
This  effective  sound  speed  is  used  directly  in  the  usual  time-dependent 
wave  Equation  (2.6). 

Since  air  behaves  very  nearly  as  an  ideal  (albeit  composite)  gas,  the 
expression  for  c  is  simply 

c2  =  y RT/H  ,  (5.3) 


where,  for  air,  y—  1.4,  R  =  8.31  J/mol*K  and  M  =  0.029  kg/mol  so  that 


c 


20.0  J  T 


(5.4) 


where  T  is  the  Kelvin  temperature.  At  this  point,  we  will  cease  to  la¬ 
bel  c  as  effective  and  understand  that  the  sound  speed  is  given  by 


c 


20.0 


+ 


U 


(5.5) 
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When  the  temperature  profile  is  such  that  convection  is  suppressed  as 
is  normally  the  case  at  night  as  a  result  of  radiative  cooling  of  the 
surface,  the  transition  layer  effects  are  primarily  surface  friction  and 
mechanically  induced  turbulence.  In  this  case,  the  PBL  is  much  thinner, 
perhaps  only  a  few  ten's  or  at  most  a  few  hundred's  of  meters.  In  this 
condition  of  stable  flow,  strong,  infrequently  turbulent  winds  aloft  are 
common  and  the  usual  positive  gradient  in  wind  speed  with  altitude  adds 
to  that  of  the  temperature  gradient  and  causes  strong  downward  refrac¬ 
tion  and  focusing  of  sound. 

In  order  to  determine  the  dominant  features  of  a  particular  PBL,  the 
stability  with  regard  to  vertical  convection  must  be  determined.  If  a 
small  parcel  of  air  is  displaced  vertically,  it  may  continue  to  move  up¬ 
ward  through  its  own  buoyancy  (in  an  unstable  layer)  or  it  may  oscillate 
about  its  rest  position  at  some  frequency  (in  a  stable  layer).  This 
frequency  is  known  as  the  Brunt -Vaisala  frequency  (for  example,  see 
Beer,  1974).  For  a  perfect  gas,  it  is 

=  g(<*  -  <**)/T  ,  (5.6) 

*■ 

where  cx  is  the  vertical  temperature  gradient,  <=<  is  the  adiabatic  gra¬ 
dient  (-0.0098  K/m) ,  and  g  is  the  acceleration  of  gravity.  If  ujg  is 
negative,  the  layer  is  unstable  and  thermal  convection  can  cause  mixing 
throughout  the  layer;  if  is  positive  then  the  flow  in  the  layer  is 
stable  and  gravity  waves  can  form. 
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In  a  convectively  unstable  PBL,  the  Richardson  number  describes  tur¬ 
bulence  maintenance  (Schlichting,  1979).  This  number  is  the  ratio  of 
buoyant  production  of  turbulence-related  kinetic  energy  to  shear-stress 

7- 

production  of  turbulent  kinetic  energy  and,  in  terms  of  c-Jg  is 

R;  =  U/<?z)2  ,  (5.7) 


where  U  is  the  horizontal  wind  speed.  If  R^  is  greater  than  about  0.25, 
turbulence  is  suppressed.  For  Rt  less  than  zero,  convective  turbulence 
dominates  and,  in  between,  shear-produced  turbulence  dominates. 

Figure  12  shows  the  meterological  parameters  that  characterize  the 
PBL  at  a  time  of  strong  focusing  (Brown,  1980  and  Thomson,  1983).  These 
parameters  were  measured  using  tethered,  instrumented  balloons  by  the 
University  of  Virginia  to  support  noise  interference  studies  for  the 
DOE/NASA  Mod-1  wind  turbine  generator  near  Boone,  North  Carolina.  The 
figure  shows  a  typical  nighttime  stable  flow.  Except  for  a  small,  ele¬ 
vated  layer  (from  75m  to  85m  altitude),  the  flow  is  stable.  In  fact,  it 
is  strongly  stable  (Rt  »  0.25)  in  the  first  50m  so  that  any  thermal 
convection  from  the  ground  would  be  strongly  suppressed.  Also,  a  strong 
vertical  wind-speed  gradient  exists  above  the  surface  layer. 

For  this  particular  PBL,  the  maximum  Brunt -Vaisala  frequency  corre¬ 
sponds  to  a  period  of  200  seconds.  Consequently,  there  will  be  no  meas¬ 
urable  coupling  between  gravity  waves  and  acoustic  waves  and  the  gravity 
terms  can  be  safely  omitted  (as  we  have  already  done)  from  the  governing 


equations . 
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Figure  13  shows  the  sound-speed  profile  for  the  meteorological  condi¬ 
tions  of  Figure  12  including  the  envelope  of  perturbations.  (This  enve¬ 
lope  is  described  in  section  5.2.)  Although  the  finite-difference  model 
can  accept  an  arbitrary  profile,  in  order  to  facilitate  comparison  to 
normal-mode  and  ray  models,  two  layers  are  used  to  fit  the  environmental 
data.  Each  layer  has  a  linear  variation  in  the  square  of  the  reciprocal 
of  the  sound  speed.  This  turns  out  to  be  an  excellent  fit  for  the  up¬ 
per,  steep-gradient  layer. 

Although  the  frequency  range  of  the  investigation  to  follow  is  so  low 
as  to  render  analysis  by  ray  theory  questionable  at  best,  some  insight 
into  the  focusing  produced  by  this  profile  can  be  gained  from  the  ray 
diagram  in  Figure  14.  The  three  limbs  of  the  caustic  (and  the  two 
cusps)  are  shown  in  Figure  15  along  with  the  observation  points  used  in 
Chapters  VI  and  VII.  Three  general  areas  are  examined  with  these 
points:  the  transition  from  the  illuminated  side  of  a  simple  caustic  to 
the  shadow  side;  the  field  between  two  closely  spaced  caustics;  and  the 
near-ground  field  before  and  near  one  of  the  cusps. 

5.2  Perturb, it  ion  Mode  1 

There  are  many  possible  processes  that  can  produce  large  scale  hori¬ 
zon:  ll  into-,  it  ios  in  the  atmospheric  boundary  layer.  One  common 
ph'-:,.  ••  •■:  •  .  . ib’.,.  conditions  is  gravity  wave  formation.  In  complex 
ter:  i  •  is  :  he  region  surrounding  the  wind  turbine  generator  at 
Bo  :.e,  ,’.;na,  mechanically  induced  turbulence  is  often  formed 
in  the  :  he  ;idge.  For  this  investigation,  we  will  choose  the 


turbulence  mechanism. 
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SOUND  SPEED  (M/S) 


Figure  13.  Sound-speed  profile  (dashed  line)  corresponding  to 

temperature  and  wind  speed  profiles  in  Figure  12.  Solid  line  shows 
profile  used  in  finite-difference  model  and  shaded  region  gives 
envelope  of  sound-speed  perturbations. 


Figure  14.  Ray  diagram  for  sound-speed  profile  in  Figure  13  (solid 
line).  Altitude  of  source  is  31m. 


Figure  15.  Location  of  ray-theory  caustics  (solid  lines)  and  groups  of 
observation  points  (crosses). 
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A  detailed  model  describing  such  turbulence  in  detail  would  be  pro¬ 
hibitively  complex  if  indeed  at  all  possible;  however,  since  we  are  not 
seeking  to  duplicate  one  particular  flow  condition  but  rather  to  de¬ 
scribe  the  general  nature  of  signal  fluctuations,  a  simple  but  not  un¬ 
realistic  model  will  suffice.  An  artificially  introduced  perturbation 
component  of  the  observed  wind  is  shown  in  Figure  16.  The  lower  part  is 
the  normal  power-law  profile  (Tennekes  and  Lumley,  1972)  typical  of  tur¬ 
bulent  boundary  flows.  This  is,  in  fact,  the  same  power- law  that  is 
used  to  extend  the  measured  profile  of  the  unperturbed  wind  to  the 
ground.  (In  other  words,  the  power-law  applies  to  the  sum  of  the  per¬ 
turbation  component  and  the  unperturbed  wind.)  This  power-law  section 
extends  to  the  top  of  the  temperature  inversion  (70m)  shown  in  Figure 
12.  Above  this  region,  there  is  no  particular  theoretical  basis  for  the 
profile  but  observations  (Thomson,  1984)  suggest  that  a  30m  layer  of 
linearly  decreasing  perturbation  speed  is  not  unreasonable.  This  per¬ 
turbation  profile  is  modulated  sinusoidally  in  range  with  a  period  of 
140m  or  twice  the  characteristic  height  and  then  added  to  the  steady- 
state  wind  profile  of  Figure  12.  In  addition,  the  perturbation  pattern 
is  advected  downwind  at  the  minimum  wind  speed  in  the  upper  part  of  the 
surface  layer  (1.5ra/s).  Since  temperature  variations  have  a  negligible 
effect  on  the  sound-speed  profile  (in  this  case),  the  steady-state  tem¬ 
perature  profile  is  used  without  perturbation. 

The  perturbed  wind  speed  is  then  given  by 


U(x,z) 


U0(z)  + 


u(z)  cos[?r(x  +  vAt)/z] 


(5.8 
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NORMALIZED  PERTURBATION  AMPLITUDE 


Figure  16. 
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where  Ua(z)  is  the  steady-state  wind  profile,  u(z)  is  the  perturbation 
amplitude  profile,  z  is  the  characteristic  height  (70m),  and  is  the 
advection  speed  (1.5m/s).  Also 

A  -  .A.'/?  A. 

|u(z/z)  ;  (0  <  z  <  z)  , 

u(z  -  z  +  30)/30  ;  (z  <  z  <  z  +  30)  ,  (5.9) 
0  ;  (z  +  30  <  z)  , 

where  u  is  the  perturbation  level  at  z.  For  this  particular  profile, 
the  wind  direction  did  not  vary  significantly  so  any  directional  correc¬ 
tions  are  ignored  and  directly  downwind  propagation  is  assumed. 

5 . 3  Source  Characteristics 

The  source  used  in  this  study  simulates  a  line  array  of  WTG's  (as 
might  be  the  case  at  a  wind  farm  rather  than  an  isolated  turbine)  by  us¬ 
ing  a  line  source.  This  is  a  crude  representation  since  the  turbines 
would  not  normally  be  rotating  in  phase;  however,  this  will  serve  to 
simulate  a  worst-case,  on-axis  beam.  In  the  finite-difference  model, 
the  line  source  is  introduced  by  using  a  point  source  in  a  twc -dimen¬ 
sional  Cartesian  coordinate  system. 

At  the  Boone  installation,  an  annoying  sound  was  produced  by  passage 
of  the  blade  through  the  turbulent  wake  of  the  supporting  tower  (Kelley, 
1981)).  To  simulate  this,  a  narrow  impulsive  source  waveform  is  used 
that  has  a  broad  spectrum  (see  Equations  (4.2)  and  (4.3)).  The  Mod-1 
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wind  turbine  generator  (WTG)  blade  rotated  at  35rpm  so  approximately  one 
pulse  per  second  was  produced.  Consequently,  there  will  be  a  rapid  drop 
in  source  level  below  1Hz. 

According  to  measurements  by  Kelley  et  al.  (1982),  the  part  of  the 
received  WTG  signal  that  was  above  the  noise  at  the  affected  dwellings 
was  bandlimited  below  25Hz.  Since  this  range  is  entirely  subaudible, 
the  annoyance  was  apparently  a  result  of  excitation  of  the  natural  vi¬ 
brational  modes  (at  approximately  5,12  and  20Hz)  of  the  human  body 
(Tempest,  1976).  The  incoming  acoustic  energy  did  not  directly  couple 
the  modes  (the  "noise"  was  not  generally  objectionable  outside  the  hous¬ 
es)  but,  instead,  it  excited  diaphragm  modes  in  the  floors  and  walls  and 
these,  in  turn,  produced  annoyances  to  the  residents. 

5.4  Configuration  for  the  PBL  Model 

Since  the  mechanics  of  the  finite-difference  method  have  been  devel¬ 
oped  in  Chapter  III,  a  particular  model  can  now  be  configured  for  the 
PBL  problem.  In  order  to  maximize  computational  efficiency,  the  range 
coordinate  is  transformed  by  gradient  scaling  (see  Appendix  B) .  In  the 
atmosphere,  as  in  the  ocean,  range  scales  are  generally  larger  than 
height  scales  but  an  approximately  square  grid  is  more  efficient  for 
calculation.  In  this  case,  a  5-to-l  range  compression  is  used  so  that 
the  resulting  grid  of  70  by  225  (2m  spacing)  represents  140m  in  altitude 
and  2250m  in  true  range. 

Because  of  this  gradient  scaling,  the  highest  sound  speed  in  the 
field  rises  to  480m/s.  In  order  to  limit  q  to  no  more  that  0.8  anywhere 
in  the  field,  the  maximum  time  step  is  0.00333s.  Notice  that,  in  the 
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unsealed  case,  this  maximum  would  be  0.00435s.  This  loss  in  time  ad¬ 
vance  per  step  is  far  outweighed  by  having  to  cover  only  one-fifth  the 
distance  downrange. 

The  source  is  located  at  31m  above  the  ground.  This  position  was  se¬ 
lected  only  to  provide  a  clean  caustic  region;  it  is  otherwise  arbi¬ 
trary.  Since  the  purpose  of  this  study  is  to  examine  only  aspects  of 
the  propagation  rather  than  simulating  the  Boone  situation  per  se,  no 
attempt  was  made  to  include  features  such  as  actual  topography  or  sur-  ** 

face  impedances . 

Source  points  traditionally  cause  problems  in  finite-difference 
schemes  (e.g.,  Alterman  and  Karal,  1968).  This  is  avoided  here  by  de¬ 
fining  a  fairly  large  source  region  (40m  in  radius  from  the  source)  and 
computing  the  source  field  analytically  at  two  consecutive  time  stops. 

These  starter  fields  are  then  used  as  initial  conditions.  A  very  accu¬ 
rate  representation  of  the  source  is  possible  if  the  peak  of  the  source 
pulse  is  put  about  halfway  out  into  this  source  region  (i.e.,  ct0  =  20 
in  this  case).  Hence,  a  start  time  of  0.059s  is  used  here. 

Dispersion  is  controlled  by  proper  selection  of  source  bandwidth  and 
grid  spacing.  These  two  parameters  determine  kh  and  this  can  be  set  ac¬ 
cording  to  the  dispersion  curve  in  Figure  6.  The  choice  of  2m  for  h  was 
made  in  order  to  keep  kh  below  3  up  to  25Hz.  Source  bandwidth  is  con¬ 
trolled  by  b  in  Equation  (4.3).  A  value  of  82.82  for  b  reduces  the 
source  amplitude  by  a  factor  of  ten  at  20Hz  and  by  a  factor  of  400  at 
30Hz . 

The  final  element  in  the  PBL  configuration  is  the  transparent-bound¬ 
ary  condition.  Since  this  investigation  is  concerned  primarily  with  re- 


fractive  focusing,  interfering  multipath  from  simple  reflection  is 
avoided  by  using  Reynolds'  (1978)  method  for  transparent  boundaries.  A 
summary  of  this  condition  is  contained  in  Appendix  A.  The  lower  limit 
on  kh  imposed  by  defects  in  the  transparent-boundary  condition  is,  ac¬ 
cording  to  Chapter  IV,  0.075,  which  is  about  2Hz  in  this  case.  Rigid 
and  pressure-release  boundaries  are,  of  course,  much  easier  to  insert. 

If  a  different  problem  is  to  be  modeled,  the  design  process  would 
proceed  as  follows: 

1.  determine  the  maximum  frequency  of  interest  and  select  b  so  that 
the  source  spectrum  drops  below  one  percent  somewhat  higher  than 
this  frequency; 

2.  based  on  the  minimum  sound  speed,  compute  the  required  grid  spac¬ 
ing  h  for  kh  less  than  T>73  at  the  maximum  frequency; 

3.  select  a  time  step  so  that  q  is  no  greater  than  0.8  at  the  maxi¬ 
mum  sound  speed; 

4.  define  a  source  region  with  a  radius  of  about  20h  and  pre-calcu- 
late  the  starter  fields  within  this  region  at  a  time  equal  to 
half  this  radius  divided  by  the  local  sound  speed; 

5.  select  appropriate  boundary  conditions. 

The  fastest  computations  are  possible  if  the  entire  solution  grid  at 
three  consecutive  time  steps  can  be  kept  in  memory,  but,  since  the  tech¬ 
nique  is  explicit,  this  is  not  a  requirement.  Pieces  of  the  grid  can  be 
swapped  in  and  out  of  active  memory  if  the  overall  grid  size  is  very 
large. 

The  model  of  atmospheric  structure  used  here  takes  about  half  a  mega¬ 
byte  of  storage  on  an  IBM  3033  so  no  memory  swapping  is  required.  One 


81 


of  the  outstanding  features  of  the  finite-difference  method  for  studies 
of  large  regions  of  the  field  is  that  the  solution  is  obtained  at  every 
grid  point,  not  just  a  few  selected  receiver  points.  A  normal-mode  so¬ 
lution  for  a  simple  two-layer  model  is  only  cheaper  than  the  finite-dif¬ 
ference  model  for  less  that  20  receiver  points. 


Chapter  VI 


STRUCTURE  OF  THE  UNPERTURBED  FIELD 

Before  the  influence  of  sound-speed  perturbations  can  be  assessed, 
the  structure  of  the  unperturbed  field  must  be  understood.  This  under¬ 
lying  structure  not  only  determines  the  nature  of  diffraction  near  a 
wave  convergence  but  also  affects  the  response  of  the  field  to  perturba¬ 
tions.  Consequently,  some  effort  is  spent  in  this  chapter  to  examine 
the  field  structure  in  the  context  of  diffractive  wave  phenomena. 

6 . 1  Overall  Field  Structure 

In  order  to  reinforce  the  distinction  between  a  strongly  diffractive 
field  and  its  asymptotic  (high-frequency)  limit  given  by  ray  theory,  let 
us  first  examine  the  overall  field  in  the  multiple  caustic  region. 

Consider  the  ray  picture  shown  in  Figure  14.  According  to  the  ray 
diagram,  the  region  of  focusing  consists  of  three  distinct  limbs  of  a 
caustic  joined  by  two  cusps.  This  structure  is  delineated  in  Figure  15. 
The  strongest  limb  is  the  ascending  limb  while  the  weakest  (by  far)  is 
the  almost-horizontal  reciprocal  limb.  (The  strengths  are  apparent  from 
the  line  density  in  the  plot  but  these  qualitative  statements  were  veri¬ 
fied  by  ray  theory  calculations  with  first-order  caustic  corrections.) 
The  almost  horizontal  limb  and  the  associated  cusp  are  called  "recipro¬ 
cal"  because  they  occur  at  about  the  reciprocal  altitude  or  the  altitude 
in  the  upper  layer  that  has  the  same  sound  speed  as  the  source  does  in 
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the  lower  layer.  A  shadow  is  formed  in  the  area  to  the  right  of  the 
ascending  limb  and  above  the  reciprocal  limb.  There  is  another  shadow 
to  the  left  and  below  the  connecting  limb  but  this  is  not  pure  shadow 
because  direct  rays  weakly  ensonify  the  region.  Also,  the  lower  cusp 
appears  to  be  particularly  intense. 

Because  the  wavelengths  considered  here  are  of  the  order  of  the  caus¬ 
tic  spacing,  conclusions  based  on  this  ray  picture  (or  on  a  ray  trace 
calculation)  can  be  quite  misleading.  To  demonstrate  this,  Figures  17, 
18,  and  19  were  constructed  from  the  results  of  the  two-layer  normal¬ 
mode  model  described  in  the  previous  chapter.  (The  finite-difference  re¬ 
sults  could  also  have  been  used.)  The  caustic  limbs  are  superimposed  on 
these  contour  plots  as  dashed  lines.  The  contour  lines  represent  lines 
of  equal  field  level  relative  to  free-field  spreading  and  the  contour 
increment  is  3dB .  In  all  cases,  the  actual  wave  field  bears  only  a  cas¬ 
ual  resemblence  to  the  ray  structure. 

The  results  at  5Hz  are  particularly  striking:  the  contoured  field  in 
Figure  17  shows  a  smooth,  broad  ridge  roughly  centered  over  the  multiple 
caustic  and  approximately  parallel  to  the  ascending  limb.  The  only  hint 
of  the  multiple  caustic  structure  is  the  slight  upturn  in  the  contour 
lines  at  the  far  right  near  the  reciprocal  limb.  There  is  no  indication 
at  all  of  the  lower  cusp  while  the  presence  of  the  upper  (reciprocal) 
cusp  merely  fills  in  the  upper  end  of  the  valley  to  the  left  of  and  par¬ 
alleling  the  focal  ridge.  The  two  major  depressions  correspond  to  the 
shadow  regions  of  the  ascending  and  connecting  limbs.  Obviously,  if  the 
caustic  structure  were  not  already  known  from  ray  tracing,  it  could  not 
be  inferred  in  any  detail  from  the  field  'ontour  plot. 


AD-A171  220 


1.1 


bi  ||M 

|so 

HIM 

m 

t  U£ 

III 22 

t  i- 

w 

■  2.0 

Ul_«. 

IS. 

1.4 

111.6 

MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BURLAU  Of  STANDARDS  1963-A 


87 


The  field  plot  at  10Hz,  shown  in  Figure  18,  delineates  the  weak  re¬ 
ciprocal  limb  more  clearly  than  at  5Hz,  but  there  is,  again,  no  indica¬ 
tion  that  there  are  two  other  strong  caustics.  Iu  fact,  at  both  5  and 
10Hz,  the  field  contour  lines  cut  across  the  connecting  limb  at  a  high 
angle  near  the  upper  cusp.  Simple  asymptotic  corrections  to  ray  theory 
would  trace  the  contours  parallel  to  the  caustic.  As  at  5Hz,  the 
"strong"  lower  cusp  does  not  show  on  the  contour  plot.  Notice  that  in 
the  upper  right  shadow  region,  the  wave  field  is  actually  highly  struc¬ 
tured  as  if  there  were  several  interfering  paths.  In  fact,  there  is 
more  structure  in  this  deep  shadow  than  in  the  illuminated  zone.  Note 
that  first-order  asymptotic  ray  theory  predicts  a  smooth,  featureless 
decay  into  the  shadow. 

At  20Hz  (Figure  19),  the  upper  cusp  and  the  wedge  formed  by  that  cusp 
and  the  connecting  and  reciprocal  limbs  are  indicated  by  the  two  ellip¬ 
tical  peaks.  Again,  the  multiple  caustic  is  represented  predominantly 
by  a  long,  single  ridge  paralleling  the  ascending  limb.  The  dropoff 
into  either  shadow  is  faster  but  the  upper  right  shadow  still  shows  much 
more  structure  than  the  illuminated  region  between  the  limbs.  Curiously 
enough,  the  lower  cusp  is  still  not  a  distinguishable  feature  of  the 
wave  field. 

By  quickly  scanning  Figures  17,  18,  and  19  in  order,  the  trend  is 
clear.  The  wave  field  is  slowly  approaching  the  general  shape  of  the 
convergence  suggested  by  the  ray  trace.  This  would  continue  until,  at 
very  high  frequency,  the  multiple  caustic  structure  would  be  clearly  de¬ 
fined  and  the  interference  pattern  in  the  shadow  would  be  generally  de¬ 
structive.  Throughout  the  frequency  range  shown  here,  however,  the  ray 
representation  is  a  very  poor  indication  of  the  actual  field. 
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In  particular,  the  wave  field  structure  is  a  generally  featureless 
enhancement  in  the  high  intensity  region  while,  in  the  shadow,  there  is 
a  distinct  interference  structure.  In  the  following  analyses,  the  caus¬ 
tics  will  be  used  as  reference  lines  to  locate  observation  points  in  the 
field,  but,  as  we  have  seen,  they  have  little  intrinsic  significance  at 
these  frequencies . 

6.2  Structure  in  the  Shadow  Zone 

Since  ray  theory  does  not  provide  much  information  about  the  field 
structure  at  low  frequency,  we  will  use  the  wave  results  given  by  both 
the  finite-difference  method  and  normal-mode  analysis.  The  shadow  zone 
is,  perhaps,  the  most  interesting  region  since  the  greatest  changes  in 
signal  characteristics  occur  there. 

6.2.1  F inite-dif ference  results 

The  finite-difference  model  produces  the  received  time  waveform  di¬ 
rectly  for  a  pulsed  source.  Figure  20  shows  the  variation  in  the  re¬ 
ceived  waveform  for  the  points  in  the  set  labelled  "shadow  zone"  in  Fig¬ 
ure  15.  The  trace  at  the  top  is  from  a  receiver  well  inside  the 
illuminated  zone  and  the  shape  is  similar  to  the  source  pulse  with  an 
extended  tail.  This  tail  is  a  characteristic  of  two-dimensional  spread¬ 
ing  and  appears  in  all  of  the  pulse  waveforms.  As  we  penetrate  the 
shadow  zone  (moving  downward  in  Figure  20),  the  pulse  flattens  consider¬ 
ably  and  broadens  slightly  as  if  the  initial  single  pulse  were  being 
split  into  two  very  closely  spaced  pulses  of  reduced  amplitude.  Another 
way  of  viewing  the  pulse  change  is  that  the  high-frequency  components 
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Figure  20.  Finite-difference  model  predictions  of  pulse  shape  at 

receiver  points  1  through  7  in  shadow  region.  Point  7  is  furthest 
into  the  classical  shadow.  Amplitude  scale  is  linear  but  otherwise 
arbitrary. 
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are  being  lost  more  quickly  deep  in  the  shadow.  As  we  will  see,  both  of 
these  interpretations  are  correct:  there  is  some  high-frequency  rejec¬ 
tion  in  the  shadow  and  there  is  a  significant  multiple  arrival  structure 
in  the  shadow  (in  spite  of  the  ray  theory  prediction  that  there  are  no 
arrivals  here  at  all). 

In  Figure  21,  the  magnitude  spectra  for  several  of  the  shadow  zone 
points  are  shown.  The  two  principal  features  are  the  increasing  amount 
of  high  frequency  rejection  with  deeper  penetration  into  the  shadow  and 
the  changing  interference  pattern.  As  we  progress  into  the  shadow,  the 
spectrum  interference  nulls  become  more  closely  spaced.  The  interfer¬ 
ence  pattern  is  quite  simple  indicating  only  two  or  three  major  inter¬ 
fering  components  and  the  null  shift  indicates  that  the  arrival  time 
difference  between  the  components  is  increasing  with  penetration  into 
the  shadow.  In  order  to  locate  the  cause  of  this  interference,  we  must 
use  another  wave  theory  model:  the  normal-mode  model. 

6.2.2  Cumulat j ve  mode -sum  analysis 

A  ray  theory  calculation  considers  the  field  to  be  the  sum  of  all 
possible  rays  connecting  the  receiver  to  the  source.  The  normal-mode 
method,  on  the  other  hand,  develops  the  field  as  the  sum  of  the  permis¬ 
sible  elementary  waves  determined  by  the  boundary  conditions.  For  a 
point  source  in  a  horizontally  stratified  medium,  the  elementary  waves 
are  cylindrical  waves  while,  in  this  analysis  for  a  line  source,  the  el¬ 
ementary  waves  are  plane  waves.  In  either  case,  if  we  cut  the  field 
with  a  vertical  plane  parallel  to  the  direction  of  propagation,  we  would 
see  that  each  mode  (elementary  wave)  would  have  its  own  wavefront  and 
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that  these  wavefronts  would  be  curved  by  the  sound-speed  gradients. 
Equivalent  rays  could  be  drawn  for  each  mode  and  they  would  be  perpen¬ 
dicular  to  the  wavefronts.  In  fact,  at  a  high  enough  frequency,  sets  of 
modes  can  be  isolated  and  summed  to  produce  a  field  pattern  that  looks 
like  a  "fuzzy"  ray  path  connecting  the  source  and  the  receiver. 

While  the  ray  plot  is  useful  for  understanding  the  structure  of  weak¬ 
ly  diffracted  (high  frequency)  fields,  the  mode-sum  plot  is  useful  for 
interpreting  the  structure  of  both  strong  and  weak  diffraction.  In  this 
plot,  the  complex  numbers  representing  each  term  in  the  mode  sum  (Equa¬ 
tion  (2.8))  are  added  and,  after  each  term  is  added,  the  new  partial  sum 
is  plotted  on  a  grid  with  the  real  and  imaginary  parts  as  the  two  axes. 
In  this  way,  the  development  or  progression  of  the  mode  sum  can  be  seen 
from  the  lowest  order  mode  through  the  highest  order  mode. 

Along  some  sections  of  this  curve,  the  modes  will  appear  to  be  com¬ 
bining  haphazardly  without  contributing  to  the  overall  movement  of  the 
endpoint  in  the  complex  plane.  Along  other  sections,  the  modes  will  be 
working  together  to  produce  generally  straight  progress  (although  often 
with  superimposed  wiggles)  in  some  direction.  In  the  first  case,  these 
modes  are  unimportant  since  they  do  not  affect  the  final  result  much. 
In  the  second  case,  the  groups  contribute  directly  to  the  field 
strength,  each  group  representing  a  specific  component  of  the  field.  At 
high  frequency,  the  groups  are  very  well  defined  and,  as  we  will  see, 
these  same  structural  features  can  be  identified  at  low  frequency  as 
well . 

Mode-sum  plots  for  four  of  the  shadow  zone  points  at  20Hz  are  shown 
in  Figure  22.  Here,  the  highest-order  modes  near  the  end  of  each  curve 
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(the  end  away  from  the  origin)  circle  around  tightly  and  do  not  contrib¬ 
ute  significantly  to  the  final  field  value.  This  is  a  good  indication 
that  enough  modes  have  been  taken  to  get  a  good  solution.  Most  of  the 
other  modes  are  adding  in  a  reasonably  orderly  manner.  Scanning  from 
the  upper  left  to  the  lower  right  (from  the  caustic  into  the  shadow),  we 
can  see  that  the  sum  seems  to  be  grouped  into  three  or  four  segments  and 
that  the  overall  curve  is  folding  up  rather  than  either  curling  up  oi 
degenerating  randomly.  Thus,  there  are,  in  the  shadow,  several  inter¬ 
fering  components  made  up  of  partial  sums  of  modes. 

If  we  compute  a  ray  equivalent  angle  based  on  the  eigenvalue  of  the 
central  mode  in  each  group  according  to  the  equation 


cos  >  (61) 

where  </>  is  the  angle  from  the  horizontal,  we  can  determine  something  of 
the  nature  of  these  groups.  In  the  upper  left  plot,  there  are  three 
distinct  groups:  one  that  goes  straight  down  from  the  origin,  a  second 
that  goes  up  and  slightly  to  the  left,  and  a  third  that  goes  further  up 
and  to  the  right.  The  second  and  third  groups  meet  at  the  real  axis. 
The  first  of  these  groups  has  an  imaginary  equivalent  angle.  This  indi¬ 
cates  an  evanescent  component  in  the  form  of  an  inhomogeneous  wave  re¬ 
sulting  from  diffractive  leakage  out  of  the  refractive  turning  point 
somewhat  below  the  receiver  altitude.  Of  course,  no  real  ray  equivalent 


exists  for  this  component. 
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The  second  and  third  groups  do  have  real  equivalent  ray  angles  and, 
thus,  represent  real  propagating  waves  that  result  from  the  diffractive 
broadening  of  the  ray  paths  forming  the  caustics  below.  A  crude  repre¬ 
sentation  of  a  low-frequency  ray  would  be  a  very  broad,  curved  swath 
and,  in  this  case,  the  paths  have  become  wide  enough  to  extend  well  into 
the  classical  shadow.  In  fact,  at  20Hz,  a  typical  value  for  the  ray 
width  parameter  given  by  Equation  (2.21)  is  100m. 

Consequently,  the  interference-like  structure  shown  in  the  shadow 
spectra  is,  in  fact,  interference  of  several  distinct  field  components 
both  propagating  and  evanescent.  The  reduction  in  field  level  into  the 
shadow  is  the  result  of  a  very  orderly  folding  of  the  mode  sum  until, 
deep  in  the  shadow,  it  has  been  broken  into  so  many  segments  that  it  ap¬ 
pears  to  be  completely  disordered. 

6.2.3  Asymptotic  analysis 

Obviously,  simple  ray  acoustics  fails  at  a  caustic;  however,  the  sin¬ 
gle-term  asymptotic  correction  given  by  Equation  (2.27)  allows  calcula¬ 
tion  of  a  value  for  the  field  at  and  near  the  caustic.  Unfortunately, 
at  the  frequency  range  of  interest  here,  this  ray  theory  correction 
gives  very  poor  results. 

Figure  23  shows  the  spectra  of  the  field  computed  by  the  finite-dif¬ 
ference  method  (solid  line)  at  points  4  and  7  of  the  shadow  zone  series. 
The  field  predicted  by  the  asymptotic  correction  to  ray  theory  is  also 
shown  as  a  dashed  line.  The  interference  structure  of  the  shadow  zone 
field  is,  of  course,  absent  in  the  asymptotic  theory. 
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This  first-order  asymptotic  correction  was  also  computed  for  the 
points  in  the  inter-caustic  zone  but  the  agreement  was  even  poorer  there 
so  that  the  results  are  not  shown.  Typical  values  for  the  characteris¬ 
tic  caustic  width  CD  of  Equation  (2.29)  range  from  50m  to  150m  at  20Hz 
(and  250ra  to  400m  at  5Hz)  so,  clearly,  the  caustics  are  much  too  close 
together  to  be  treated  as  distinct.  Consequently,  the  first-order  as¬ 
ymptotic  correction  is  of  little  value.  We  will  consider  the  second-or¬ 
der  correction  in  the  following  section. 

6.3  Structure  of  the  Multiple  Caustic 
6.3.1  Mode -sum  analysis 

As  in  the  shadow  zone,  we  will  again  resort  to  the  mode-sum  plots  to 
interpret  the  field  structure.  Since  there  is  very  little  change  in  the 
spectrum  or  waveshape  of  the  received  pulses  over  the  inter-caustic  zone 
(and  the  cusp  zone),  the  spectra  and  time  series  plots  will  not  be 
shown.  Figure  24  shows  the  mode-sum  plots  for  four  of  the  points  in  the 
inter-caustic  region  of  Figure  15.  The  lower  left  plot  is  at  the  point 
of  maximum  field  level  in  this  zone. 

One  feature,  identical  (except  for  angular  orientation)  in  all  four 
of  the  plots,  is  that  of  the  first  two  straight  line  segments  from  the 
origin.  These  two  modes  give  the  direct,  low-angle  path  into  this  zone 
that  is  also  shown  on  the  ray  plot.  The  other  modes  are  associated  pre¬ 
dominantly  with  the  higher  angle  converging  rays. 

The  lower  left  plot  at  the  field  maximum  shows  that  all  of  the  modes 


form  one  group  (beyond  the  direct  path  segment  that  happens  to  be  in  the 
same  direction)  working  together  to  build  up  the  field  level.  Instead 
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of  breaking  up  into  a  few  discrete  groups  as  the  shadow  zone  modes  did, 
this  single  mode  group  smoothly  curls  up  on  either  side  of  the  field 
maximum.  There  is  no  separable  two-caustic  structure  at  this  frequency 
(20Hz) .  This  is  in  accord  with  the  overall  field  levels  shown  in  Figure 
19.  Notice  at  point  1  (upper  left  plot  of  Figure  24)  the  mode  sum  has 
curled  up  to  form  a  complete  circle  without  showing  any  other  distinct 
grouping. 

Hence,  the  field  in  the  inter-caustic  zone  actually  only  has  one 
structural  feature  (aside  from  the  trivial  direct  path)  and  this  is  rep¬ 
resentative  of  the  general  convergence  and  constructive  interference  in 
this  area.  There  is  no  indication,  even  in  the  mode  sum,  that  there  are 
two  caustics.  Consequently,  use  of  the  extended  yl  theory  to  predict 
fluctuations  in  the  inter-caustic  zone  probably  has  little  meaning  at 
these  frequencies  since  it  takes  both  caustics  into  account. 

6.3.2  Asymptotic  analysis 

As  we  have  seen,  the  two  caustics  bordering  the  inter-caustic  zone 
cannot  be  treated  as  distinct  elements.  Instead,  we  can  compute  the 
field  based  on  the  second-order  asymptotic  correction  to  ray  theory  giv¬ 
en  in  Equation  (2.31).  Figure  25  gives  a  contour  plot  of  the  field  lev¬ 
el  (relative  to  free-field  spreading)  predicted  by  this  second-order 
correction.  To  construct  this  figure,  the  weak  reciprocal  limb  was  ig¬ 
nored  and  the  range  Equation  (2.30)  was  fit  to  the  ascending  and  con¬ 
necting  limbs.  These  two  limbs  did,  however,  have  quite  different  sec¬ 
ond  derivatives  at  the  caustics  so  a  symmetric  cubic  curve  is  not  a  good 
representation.  This  asymmetry  was  reduced  greatly  by  bisecting  the 


Figure  25.  Prediction  of  field  level  at  20Hz  by  second-order  asymptotic 
correction  to  ray  theory.  Compare  to  Figure  19. 
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wedge  formed  by  these  two  caustics  and  the  lower  cusp,  and  taking  a  new 
"range"  coordinate  perpendicular  to  this  bisector.  In  this  way,  the 
transformed  range  coordinate  cuts  both  caustics  at  an  angle  approaching 
90  degrees.  The  second  derivatives  of  this  transformed  range  became 
much  more  nearly  equal.  Once  the  coordinate  transformation  was  done,  X 
and  Y  could  be  computed  by  Equations  (2.33)  for  a  grid  superimposed  on 
the  caustic  system.  Then  the  field  was  calculated  from  Equation  (2.31) 
as  discussed  in  Chapter  II. 

Although  the  frequency  is  still  too  low  to  allow  detailed  prediction 
of  absolute  field  levels,  the  desired  structure  is  clearly  shown.  There 
is  a  smooth  ridge  centered  on  the  two  caustics  with  depressions  on  ei¬ 
ther  side.  In  addition,  beyond  these  depressions,  the  field  temporarily 
increases  again  thus  showing  the  structure  that  is  observed  in  the  shad¬ 
ow  by  the  mode  and  finite-difference  methods.  This  structure  is  com¬ 
pletely  absent  from  the  first-order  asymptotic  calculations. 

As  the  frequency  is  increased,  the  major  depressions  would  slide  in 
toward  the  caustics  and  down  toward  the  lower  cusp.  By  scanning  Figures 
17,  18,  and  19,  we  can  see  that  this  trend  is  correct  but,  again,  the 
frequency  is  not  high  enough  for  the  correspondence  to  be  good  in  de¬ 
tail. 

6.4  Structure  near  the  Cusp 

Because  this  cusp  does  not  appear  as  a  distinct  feature  on  the  over¬ 
all  field  plot  in  Figure  19,  we  would  not  expect  the  structure  to  be 
very  different  from  that  in  the  inter-caustic  zone. 
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The  mode-sum  plots  for  four  points  of  the  cusp  series  (see  Figure  15) 
are  shown  in  Figure  26.  As  in  the  inter-caustic  zone,  there  are  two 
distinct  features:  one  mode  representing  the  low-angle  direct  path  and  a 
large  group  of  modes  that  represents  the  convergence.  The  direct  path 
mode  here  generally  opposes  the  other  group  but,  otherwise,  these  plots 
are  very  similar  to  those  in  Figure  24 . 

We  would  expect  that  a  frequency  high  enough  to  display  the  multiple 
caustic  structure  in  the  inter-caustic  region  would  also  resolve  the 
cusp.  By  means  of  the  X  parameter  of  the  second-order  correction,  we 
can  determine  what  the  lower  frequency  limit  for  resolution  of  the  two 
caustics  is.  A  value  of  -4  for  X  is  required  in  order  for  the  integral 
of  Equation  (2.31)  to  show  two  distinct  ridges,  one  at  each  caustic.  At 
20Hz,  X  is  about  -2  at  50m  altitude.  Since  the  horizontal  wave  number 
is  proportional  to  frequency,  each  derivative  of  range  with  respect  to 
wave  number  gives  an  inverse  proportionality  to  frequency.  Thus,  from 
the  first  of  Equations  (2.33),  we  can  see  that  X  is  proportional  to  the 
square  root  of  frequency.  Consequently,  the  two  caustics  and  the  cusp 
should  appear  as  distinct  features  at  80Hz  (four  times  20Hz)  or  greater. 


Chapter  VII 
SIGNAL  FLUCTUATIONS 

Once  time-varying  perturbations  in  sound  speed  are  introduced,  the 
signal  levels  throughout  the  field  fluctuate.  These  fluctuations  are 
easy  to  simulate  with  the  range-dependent  finite-difference  model,  but 
this  simulation,  by  itself,  does  not  shed  much  light  on  the  underlying 
causes  of  the  fluctuations. 

In  Chapter  III,  we  considered  a  ray  interpretation  of  the  fluctuation 
mechanism  through  the  and  gT  parameters.  By  following  rays  through  a 
perturbed  medium,  three  effects  can  be  isolated:  First,  the  ray  path 

itself  changes  from  the  unperturbed  path  and  there  is  a  concommittant 
change  in  spreading  loss.  Second,  an  additional  path  may  be  possible, 
thereby  introducing  the  potential  for  multipath  interference.  These  two 
effects  are  considered  simultaneously  by  the  diffraction  parameter 
and,  in  cases  where,  under  the  influence  of  perturbations,  the  original 
path  splits  into  two  paths,  neither  of  which  resembles  the  original 
path,  the  /[.  parameter  can  miss  important  contributions  to  the  fluctua¬ 
tion.  The  third  effect  that  is  described  by  ray  theory  is  the  phase 
change  resulting  from  differences  in  sound  speed  along  the  path.  This 
effect  is  given  by  the  strength  parameter  <p. 

When  the  unperturbed  field  intensity  can  be  adequately  modeled  by  ray 
theory,  these  ray-based  concepts  of  fluctuation  have  some  utility.  How¬ 
ever,  by  probing  a  strongly  diffractive  region  of  focusing  by  wave  meth- 
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ods  in  Chapter  VI,  we  have  seen  that  ray  theory  predictions  can  recreate 
the  general  features  of  the  field  only  after  a  rather  elaborate  asymp¬ 
totic  correction.  Since  the  fluctuation  parameters  A.  and  are  based 
on  simple  ray  theory,  there  is  no  reason  to  expect  that  these  parameters 
will  be  useful  for  estimating  the  type  and  magnitude  of  signal  fluctua¬ 
tions  in  the  present  example. 

In  order  to  characterize  fluctuations  in  wave  terms  then,  we  will 
consider  the  normal-mode  decomposition  of  the  field.  Ihis  approach  is 
still  deficient  in  that,  in  the  case  we  are  considering,  mode  coupling 
is  undoubted^  significant  and  we  will  ignore  this  coupling  completely. 
In  spite  of  this  omission,  some  useful  results  can  be  developed.  In 
this  chapter,  the  fluctuation  predictions  of  the  finite-difference  model 
will  be  briefly  summarized  and  then  a  substantial  effort  will  be  devoted 
to  the  mode  interpretation  of  these  fluctuations. 

7 . 1  Calculated  Fluctuations 

As  described  in  Chapter  V,  the  perturbations  are  introduced  as  a  hor¬ 
izontally  cyclic  pattern  of  sound-speed  variation  that  drifts  through 
the  field  at  1.5m/s  as  a  crude  representation  of  the  advection  of  large- 
scale  turbulence.  The  pattern's  period  is  140m  so  the  fluctuations 
should  completely  cycle  about  every  95s.  As  a  compromise  between  the 
expense  of  running  the  finite-difference  model  and  the  desirability  of 
frequent  sampling,  sixteen  runs  were  made  covering  one  complete  cycle  of 
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7.1.1  Time  history  of  received  levels 

Figures  27  and  28  show  the  calculated  level  fluctuations  at  2,  5,  10 
and  20Hz  for  two  of  the  points  in  the  inter-caustic  zone  (points  1  and 
5).  Here,  as  before,  the  relative  level  scale  indicates  the  gain  over 
free-field  spreading  in  a  homogeneous  medium.  As  we  would  expect,  the 
higher  frequencies  show  greater  fluctuation  and  also  greater  departure 
from  free-field  spreading.  Notice  also  that  the  20Hz  curve  shows  more 
variability  at  point  1  than  at  point  5,  point  5  being  at  a  point  of 
greater  field  enhancement.  This  is  a  consistent  feature  of  these  simu¬ 
lations:  in  regions  of  strong  focusing,  the  signal  variability  (at  the 
higher  frequencies)  is  generally  lower  than  in  other  regions.  This  be¬ 
havior  is  modeled  in  section  7.3. 

The  behavior  of  the  signal  level  near  the  cusp  as  a  function  of  time 
is  not  significantly  different  from  the  behavior  in  the  inter-caustic 
zone  but  there  are  some  noteable  differences  in  the  shadow  zone.  The 
signal-level  time  histories  for  points  2  and  5  in  the  shadow  region  are 
shown  in  Figures  29  and  30,  respectively.  The  progression  in  relative 
gain  with  frequency  is  not  as  smooth  as  in  the  inter-caustic  zone  but 
this  can  also  be  seen  in  the  unperturbed  spectra.  As  before,  the  higher 
frequencies  have  greater  level  changes  and  the  greatest  changes  seen  in 
the  20Hz  curve  in  Figure  30  are  accompanied  by  a  marked  asymmetry.  This 
will  be  considered  in  more  detail  later  but,  for  now,  it  is  enough  to 
say  that  this  is  an  indication  of  alternately  constructive  and  destruc¬ 
tive  interference  between  two  nearly  equal  components. 


Figure  27.  Signal  level  fluctuations  at  point  1  in  the  inter-c3ust  ic 

zone.  Frequencies  shown  are  2  ( - -  — ),  5  ( - ),  10  ( - 

and  20Hz  ( - ). 
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7.1.2  Normalized  fluctuation  magnitude 

In  order  to  compare  the  fluctuations  at  different  receiver  locations 
and  frequencies,  the  normalized  fluctuation  magnitude  V  given  by  Equa¬ 
tion  (3.12)  will  be  used.  This  quantity  is  actually  the  standard  devia¬ 
tion  of  the  received  amplitude  divided  by  the  average  received  ampli¬ 
tude.  The  fluctuation  magnitude  is  a  convenient  quantity  because  it  is 
a  measure  of  signal  variability  that  is  independent  of  the  mean  signal 
level.  Also,  if  the  signal  is  determined  by  the  interference  of  many 
components  of  random  amplitude  and  phase,  V  approaches  an  upper  limit  of 
0.52  (Skudrzyk,  1957).  This  is  the  saturation  limit.  (If,  however,  the 
signal  results  from  interference  of  only  two  components,  the  limiting 
value  of  V  can  be  somewhat  less.  This  will  be  demonstrated  in  section 
7.3.) 

The  values  for  fluctuation  magnitude  for  receivers  in  all  three  of 
the  observation  areas  and  at  2,  5,  10  and  20Hz  are  shown  in  Figure  31. 
While  the  20Hz  values  are  consistently  higher  than  all  other  values, 
there  is  no  other  definite  trend  with  respect  to  frequency.  None  of  the 
values  exceed  the  saturation  limit  of  0.52  although  the  20Hz  curves  ap¬ 
proach  this  value.  As  we  will  see  in  the  next  section,  the  20Hz  compo¬ 
nent  of  the  signal  is  just  beyond  saturation  at  the  ranges  considered 


here . 


Figure  31.  Fluctuation  magnitude  computed  by  the  finite-difference 
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7 . 2  Mode  Analysis  of  Fluctuat ions 

From  Figure  31  it  is  clear  that  the  normalized  fluctuation  magnitude 
does  not  follow  a  simple  power  law  in  range  as  is  predicted  by  the  homo¬ 
geneous,  isotropic  scattering  model.  The  field  in  this  region  is  a  com¬ 
plicated  combination  of  coherent  components  and  any  attempt  to  analyze 
the  nature  of  fluctuations  must  account  for  this  structure.  The  analy¬ 
sis  to  follow  is  based  on  normal-mode  theory  and,  in  this  context,  the 
following  points  are  important: 

1.  the  degree  of  influence  of  the  perturbations  on  a  given  mode; 

2.  the  relative  contribution  of  that  mode  to  the  received  field  (or, 
in  other  words,  the  degree  to  which  that  mode  is  excited); 

3.  the  structure  of  the  field  at  a  particular  receiver;  and, 

4.  the  influence  of  the  range  distribution  of  the  perturbations. 

The  first  three  points  can  easily  be  described  by  mode  theory  but  the 
fourth  is  more  difficult  unless  the  modes  are  allowed  to  couple,  in 
which  case  the  analysis  becomes  more  time-consuming  than  the  complete 
finite-difference  simulation.  This  defeats  the  purpose  of  the  modal 
analysis,  which  is  to  provide  a  simple  model  for  the  fluctuations  with¬ 
out  solving  the  complete  field  problem.  Consequently,  mode  coupling 
will  be  neglected.  In  spite  of  this  omission,  a  fairly  good  estimate  of 
the  fluctuation  magnitude  can  be  made  and  some  of  the  important  fluctua¬ 


tion  mechanisms  can  be  identified. 
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7.2.1  Equivalent  fluctuation  magnitude 

By  using  the  theory  developed  in  section  3.2.3,  we  can  estimate  the 
eigenvalue  shift  introduced  by  some  perturbation  in  sound  speed.  To  do 
this,  we  must  first  fit  the  perturbation  function  given,  in  this  case, 
by  Equation  (5.9)  with  functions  that  allow  evaluation  of  the  eigenvalue 
shift  expression,  Equation  (3.29). 

There  is,  of  course,  no  point  in  resolving  the  perturbation  profile 
in  much  finer  detail  than  a  wavelength  at  the  highest  frequency  of  in¬ 
terest.  Surprisingly  good  results  are  obtained  by  fitting  Equation 
(5.9)  with  only  two  segments:  one  segment  linear  in  altitude  above  an¬ 
other  constant  segment.  By  means  of  Equation  (3.34),  these  segments  can 
be  written  as 


£c/c3  =  15  [l  -  (z  -  58)/42]  /3383  ,  (7.1) 


for  z  between  58m  and  100m,  and 

£c/cJ  =  15/3383  ,  (7.2) 

for  z  between  0m  and  58m.  Otherwise,  the  sound-speed  perturbation  is 
zero. 

By  using  these  expressions  in  Equation  (3.29),  we  can  calculate  the 
maximum  positive  and  negative  eigenvalue  shifts.  These  shifted  eigenva¬ 
lues  are  then  used  in  the  exponent  of  the  mode-sum  Equation  (2.8)  to 
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compute  the  upper  and  lower  limits  for  the  fluctuation  pressure  ampli¬ 
tude.  These  two  values  and  the  unperturbed  pressure  amplitude  are  used 
to  compute  V  by  means  of  Equation  (3.12).  Since  the  amplitude  fluctua¬ 
tions  are  usually  roughly  sinusoidal,  this  over-estimates  V  so  the  val¬ 
ues  are  corrected  by  assuming  that*  the  extreme  values  of  pressure  ampli¬ 
tude  are  the  upper  and  lower  peaks  of  a  sinusoidal  curve. 

The  results  of  these  calculations  for  5  and  20Hz  are  shown  in  Figure 
32.  The  solid  lines  reproduce  the  finite-difference  calculations,  while 
the  dashed  lines  give  the  calculations  by  eigenvalue  shift.  The  results 
at  2  and  10Hz  are  generally  better  in  the  shadow  region  and  worse  in  the 
inter-caustic  zone. 

There  are  three  obvious  flaws  in  the  eigenvalue  shift  calculation: 
First,  there  is  no  attempt  to  consider  the  actual  range  distribution  of 
the  perturbations.  This  does  not  seem  to  be  a  very  serious  problem 
since  the  correspondence  between  the  eigenvalue  shift  calculation  and 
the  finite-difference  calculation  is  quite  good  over  the  entire  inter¬ 
caustic  zone .  '■ 

Second,  because  the  shadow  region  has  significant  interference  struc¬ 
ture,  the  actual  location  of  peaks  and  nulls  can  be  sensitive  to  small 
changes  in  the  profile  or  boundary  conditions.  Recall  that  the  normal¬ 
mode  model  uses  an  extended  sound-speed  profile  instead  of  explicit 
boundary  conditions;  hence,  the  structure  in  the  shadow  zone  is  repro¬ 
duced  accurately  in  form  but  not  necessarily  in  precise  position.  As  a 
result,  a  point  for  which  the  finite-difference  model  predicts  an  inter¬ 
ference  peak  may,  according  to  the  normal-mode  model,  be  near  a  null. 
This  displacement  in  structural  details  can  be  seen  by  comparing  the 
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20Hz  values  of  the  finite-difference  spectra  in  Figure  21  with  the  nor¬ 
mal-mode  field  plot  in  Figure  19.  This  is  far  more  likely  to  be  a  prob¬ 
lem  at  the  higher  frequencies  and  we  do  see  considerable  disagreement  in 
curve  shape  at  20Hz  in  Figure  32  for  the  shadow  region  points.  The 
curve  shape  at  5Hz  in  the  shadow  region  is  reproduced  accurately  but 
with  a  shift  in  level.  The  cause  for  this  is  not  clear  but  it  probably 
is  not  related  to  the  structure  problem  at  20Hz.  The  2  and  10Hz  curves 
(not  shown)  in  the  inter-caustic  zone  also  reproduce  the  required  shape 
accurately  but  with  a  significant  level  shift. 

The  final  problem  with  the  eigenvalue  shift  calculation  results  from 
the  assumption  that  the  maximum  eigenvalue  shift  will  produce  the  maxi¬ 
mum  level  fluctuation.  This  is  reasonable  only  below  saturation.  Above 
some  frequency  (or  beyond  some  range),  the  mode-sum  pattern  will  be  dis¬ 
turbed  so  much  that  the  magnitude  of  the  sum  will  appear  to  vary  random¬ 
ly  rather  than  continue  to  increase  or  decrease.  As  we  will  see  short¬ 
ly,  the  points  considered  here  do  show  saturation  at  20Hz  and  the 
eigenvalue  shift  calculations  could  be  affected  by  being  too  low.  The 
disagreement  in  the  vicinity  of  the  cusp  is  more  likely  a  result  of  the 
proximity  of  the  boundary  since  the  boundary  is  modeled  differently  by 
each  technique. 

7.2.2  Modal  ana lys is 

We  now  have  two  methods  for  calculating  the  fluctuation  magnitude  but 
neither  one  gives  much  physical  insight  into  the  fluctuation  process. 
For  this,  we  must  examine  the  components  of  the  mode  sum  in  the  context 
of  the  first  three  points  listed  at  the  beginning  of  section  7.2,  name- 
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ly:  the  perturbation- induced  change  in  a  given  mode;  that  mode's  rela¬ 
tive  contribution  to  the  sum;  and  the  overall  field  structure. 

By  computing  the  eigenvalue  shift,  we  have  already  found  the  change 
in  the  individual  modes.  Before  examining  this  more  closely,  consider 
another  approach  to  finding  the  degree  of  influence  of  the  perturbations 
on  the  modes:  expand  the  perturbation  function  Jc/c*  as  a  sum  of  the 
unperturbed  eigenfunctions.  The  series  coefficient  of  any  eigenfunction 
is  then  the  relative  excitation  of  that  eigenfunction  by  the  perturba¬ 
tion  function.  Since  the  eigenfunctions  are  orthogonal,  the  coeffi¬ 
cients  for  the  series  expansion  are  easily  computed  and  the  result,  to  a 
constant  factor,  is  identical  to  Equation  (3.29),  the  equation  for  the 
eigenvalue  shift.  (We  have  also  used  Equation  (3.34),  valid  for  small 
perturbations.)  Therefore,  the  series  coefficients  are  simply  the  ei¬ 
genvalue  shifts. 

Now,  the  difference  between  the  eigenvalues  before  and  after  pertur¬ 
bation  does  not  directly  indicate  the  change  in  the  appropriate  term  of 
the  mode  sum.  Since  we  are  considering  only  the  phase  changes  in  the 
mode  sum  to  be  important,  we  would  like  to  know  the  difference  between 
the  old  phase  and  the  perturbed  phase.  This  depends  not  only  on  the  ei¬ 
genvalue  shift  but  also  on  the  horizontal  range.  Furthermore,  this 
phase  difference  can  be  greater  than  a  full  cycle  but  we  are  only  inter¬ 
ested  in  the  net  change  in  the  complex  terra  in  the  mode  sum.  Conse¬ 
quently,  let  us  define  a  perturbation  weight  H  such  that 


sin[(yV  ■  >x/2] 


(7.3) 
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When  H  is  zero,  the  eigenvalue  shift  is  such  that  there  is  no  net  change 
in  the  mode  phase;  when  H  is  unity,  the  phase  has  shifted  by  180  de¬ 
grees.  In  other  words,  the  larger  H  is  for  a  given  mode,  the  more  that 
mode  is  changed  by  the  perturbations. 

A  sample  calculation  of  this  perturbation  weight  is  shown  in  Figure 
33  for  points  1  and  5  in  the  inter-caustic  zone  at  5Hz.  The  results  for 
a  similar  calculation  at  20Hz  are  shown  in  Figure  34.  For  the  particular 
shape  of  perturbation  function  used  in  this  study,  the  lowest  mode  is 
always  the  most  affected  mode,  although  the  phase  change  may  be  greater 
than  180  degrees,  in  which  case  the  perturbation  weight  drops  below  uni¬ 
ty.  This  is  the  case  for  the  solid  curve  in  Figure  34.  If  the  frequen¬ 
cy  is  high  enough  (or  the  range  large  enough),  the  behavior  of  H  for  the 
lower  modes  can  appear  to  be  quite  erratic  since  the  total  phase  differ¬ 
ence  can  be  several  complete  cycles.  This  is  illustrated  by  Figure  35 
in  which  the  calculation  of  the  previous  figures  is  repeated  at  50Hz. 
Here,  the  first  six  modes  show  large  phase  changes  over  the  range  inter¬ 
val  from  point  1  to  point  5. 

By  scanning  Figures  33  through  35,  we  can  see  two  distinct  kinds  of 
perturbation  weight  changes  with  range.  The  first  kind  of  change  is 
seen  in  the  low-order  modes,  which  show  the  most  change  with  perturba¬ 
tion  at  any  given  range  and  also  show  the  most  change  as  a  function  of 
range.  (This  range  dependence  is  only  clear  at  the  higher  frequencies.) 
The  second  kind  of  change  is  the  relatively  small  but  consistent  change 
between  points  1  and  5  (i.e.,  between  the  solid  and  the  dashed  curves) 
for  the  modes  beyond  the  low-order  group.  This  change  increases  with 
frequency  and  slowly  decreases  with  mode  number.  Since  the  phase  change 
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Figure  34.  Perturbation  weight  function  for  the  normal-mode  solutions 
at  points  1  ( - )  and  5  ( - )  in  the  inter-caustic  zone  (20Hz)  . 


Figure  35.  Perturbation  weight  function  for  the  normal-mode  solutions 
at  points  1  ( - )  and  5  ( - )  in  the  inter-caustic  zone  (50Hz). 
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is  relatively  constant  with  mode  number,  the  part  of  the  mode  sum  corre¬ 
sponding  to  higher  order  modes  retains  its  shape  and  rotates  as  a  whole 
by  that  small  phase  change. 

At  the  lower  frequencies,  the  perturbations  affect  only  the  low-order 
modes  to  any  significant  degree.  If  these  modes  are  important  to  the 
sum  (we  will  see  later  how  to  determine  this),  then  the  fluctuations  may 
be  large.  At  the  higher  frequencies,  the  perturbation  effects  on  the 
low-order  modes  increase  until  saturation  is  reached.  The  onset  of  sat¬ 
uration  is  indicated  when  the  phase  shift  of  the  most-affected  mode  goes 
beyond  180  degrees.  From  Figure  34  and  35,  we  can  see  that  the  fluctua¬ 
tions  are  just  saturated  at  20Hz  and  well  into  the  saturation  region  at 
50Hz . 

Now  that  we  have  seen  how  the  perturbations  affect  individual  modes, 
we  will  consider  how  the  individual  mode  changes  influence  the  total 
field  through  the  mode  sum.  Not  all  modes  contribute  equally  to  the 
field;  their  contribution  depends  on  how  strongly  they  are  excited  and 
this,  in  turn,  depends  on  the  location  of  the  source  and  receiver  with 
respect  to  the  nodes  or  antinodes  of  the  mode  eigenfunctions.  If  we 
call  the  amplitude  of  any  one  mode  in  the  mode  sum  given  by  Equation 
(2.8)  its  excitation  and  plot  this  excitation  as  a  function  of  mode  num¬ 
ber,  we  obtain  a  curve  like  that  in  Figure  36.  Since,  for  uncoupled 
modes,  the  excitation  is  not  a  function  of  range,  only  of  source  and  re¬ 
ceiver  altitude,  this  one  curve  applies  to  all  of  the  inter-caustic  zone 
points  at  5Hz.  The  same  calculation  at  20Hz  is  shown  in  Figure  37.  We 
can  see  that,  at  both  frequencies,  the  perturbation  changes  of  the  low- 
order  modes  will  be  important  since  the  low-order  modes  are  the  most 
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The  results  in  the  shadow  region  are  somewhat  different.  Figure  38 
shows  that  the  first  few  modes  for  a  point  near  the  focus  (point  2)  are 
strongly  excited  while  deeper  in  the  shadow  (at  point  5),  the  higher 
modes  are  stronger.  This  feature  is  consistent  at  all  frequencies  and, 
therefore,  allows  us  to  determine  how  much  of  an  influence  changes  in 
the  low-order  modes  have  on  the  fluctuation  magnitude.  Referring 
back  to  Figure  32,  we  can  see  that,  at  5Hz  in  the  shadow  region  (and, 
also,  for  the  10Hz  curve  not  shown),  the  trend  is  toward  lower  fluctua¬ 
tion  magnitude  deeper  in  the  shadow.  Hence,  here,  the  varying  excita¬ 
tion  of  the  low-order  modes  is  controlling  the  fluctuation  magnitude. 
At  20Hz,  there  is  no  similar  trend. 

Finally,  the  degree  of  signal  fluctuation  also  depends  on  the  struc¬ 
ture  of  the  field.  As  demonstrated  by  the  mode-sum  curves  in  Chapter 
VI,  the  structure  in  the  shadow  region  is  substantially  different  from 
that  in  the  inter-caustic  zone.  In  between  the  caustics,  where  the 
field  is  strong,  the  modes  add  in  a  very  orderly  manner  (for  example, 
see  Figure  24,  the  lower  left  plot).  Here,  we  would  expect  the  fluctua¬ 
tions  to  be  more  predictable  and,  in  fact,  the  agreement  between  the  fi¬ 
nite-difference  calculation  of  V  and  the  eigenvalue  shift  calculation  is 
quite  good  in  this  region. 

7 . 3  Two -component  Model 

Particularly  in  the  region  of  strong  enhancement  between  the  caus¬ 
tics,  we  have  been  able  to  successfully  predict  the  fluctuation  magni¬ 
tude  by  computing  the  eigenvalue  shift  introduced  by  sound-speed  pertur¬ 
bations  and  then  summing  the  phase-shifted  modes.  While  this  method 
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avoids  solving  the  range-dependent  problem  directly,  its  principal  value 
rests  in  the  insight  gained  into  the  fluctuation  process.  In  this  sec¬ 
tion,  we  will  construct  a  very  simple  fluctuation  model  based  on  just 
one  of  the  properties  observed  in  the  eigenvalue-shift  work. 

In  Chapter  V,  we  saw  that  the  field  structure  in  the  shadow  region 
comprised  two  or  three  dominant  components,  which  produced  a  rather  sim¬ 
ple  interference  pattern.  Furthermore,  the  fluctuating  signal  level  at 
one  of  the  points  in  the  shadow  (see  Figure  30,  20Hz  curve)  appeared  to 
be  the  result  of  two  nearly  equal  components  with  variable  relative 
phase.  Finally,  the  eigenvalue  shift  always  affected  a  few  of  the  low¬ 
est  order  modes  significantly  more  than  the  other  modes  (see  Figures  33 
and  34).  These  observations  suggest  that  a  model  consisting  of  the  vec¬ 
tor  sum  of  two,  unequal  amplitude  components  may  be  able  to  describe  the 
signal  level  fluctuations  if  the  relative  phase  of  these  components  is 
varied  according  to  the  eigenvalue  shift  calculations. 

We  have,  then,  a  sum  of  two  vectors  of  magnitude  P„  and  crPB  with 
relative  phase  ^ .  The  resultant  magnitude  is 

P  =  Pc  yi  +  cr*  +  2  a  cos  ,  (7.4) 


and  the  average  resultant  for  ^  varying  over  some  range  from  </  to 
is 


P 


(7.5) 
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or 


P 


2(1  +tf)  E(,4/2,r) 


E(^/2,r) 


(7.6) 


where  E  is  the  incomplete  elliptic  integral  of  the  second  kind  and 


r 


+  a) 


(7.7) 


The  resulting  fluctuation  magnitude  can  be  computed  by  integrating  to 
find  the  mean-square  amplitude  difference  (see  Equation  (3.12)) 


1  + 


1  +  a- 


1  +  2  c?  (s in  jZ^  -  sin  ^  )/  (  <£z  -  )J 


(7.8) 


(P0/P)2 


7.3.1  Pred ict ions  at  saturat ion 

The  maximum  fluctuation  range  results  when  varies  from  -  r  to  +7r 
(or  more,  of  course);  this  is  the  saturation  limit  for  this  model.  Be¬ 
cause  of  symmetry,  we  can  calculate  the  saturation  limit  by  setting  ^ 
equal  to  zero  and  equal  to  'Tr'.  Abramowitz  and  Stegun  (1965)  give  a 
convenient  polynomial  approximation  for  E(7>V 2,r)  so  that  P  may  be  cal- 
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culated  easily.  If  the  two  components  are  equal  (that  is,  cj"  =  1),  then 
the  fluctuation  magnitude  V  is  equal  to  0.48.  This  is  an  absolute  limit 
on  V  for  this  simple  model  and  it  agrees  quite  well  with  the  usual  satu¬ 
ration  limit  of  0.52  (Skudrzyk, 1957) .  Figure  39  gives  the  relationship 
between  V  and  o'.  As  we  can  see,  if  the  two  dominant  field  components 
are  unequal  in  amplitude,  the  fluctuation  magnitude  can  be  well  below 
0.52  even  at  saturation. 

As  a  check  of  the  validity  of  the  two-component  model,  let  us  take 
the  minimum  and  maximum  amplitudes  of  the  fluctuation  time  history 
curves  for  20Hz  and  compute  the  fluctuation  magnitude  using  Equation 
(7.8).  If  we  take,  in  particular,  the  curves  corresponding  to  points  2 
and  5  in  the  shadow  region  (Figures  29  and  30),  we  find  that  the  minimum 
and  maximum  levels  at  point  2  are  -4.1dB  and  -l.OdB  respectively  while 
those  at  point  5  are  -25.9dB  and  -8.9dB.  Consequently,  o'  is  0.18  for 
point  2  and  0.75  for  point  5.  The  corresponding  fluctuation  magnitudes 
computed  by  Equation  (7.8)  are  then  0.12  and  0.43,  respectively.  These 
values  are  essentially  identical  to  those  computed  by  the  finite-differ¬ 
ence  model  and  shown  in  Figure  32.  While  this  demonstrates  the  consis¬ 
tency  of  the  two-component  method,  the  fluctuation  time  history  minima 
and  maxima  are  not  easy  to  calculate  and  this  is  not  a  practical  way  to 
compute  V. 

To  apply  the  two-component  model  to  the  problem  considered  here,  a 
practical  rule  for  selecting  the  two  components  must  be  established. 
The  primary  indicator  for  separating  the  two  components  is  the  relative 
shift  in  eigenvalue  of  the  modes  as  given  by  the  perturbation  weight 
function.  (See,  for  example,  Figure  34.)  For  the  particular  perturba- 
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tion  function  used  in  this  analysis,  the  shift  is  always  greatest  for 
the  low-order  modes.  The  eigenvalues  of  the  modes  near  this  peak  in  the 
perturbation  weight  curve  have  approximately  the  same  shift;  these 
shifts  are  correlated  with  respect  to  the  perturbations.  This  can  be 
seen  in  a  qualitative  way  by  examining  the  three  lowest  modes  in  Figure 
24.  Although  they  rotate  in  the  plane  from  diagram  to  diagram,  they  ro¬ 
tate  as  a  group  with  no  significant  shift  relative  to  one  another. 

At  both  5  and  20Hz,  the  first  component  would  consist  of  several  of 
the  lowest  modes;  however,  not  all  of  these  modes  are  significant.  The 
mode  amplitude  or  excitation  curve  (see  Figure  37)  drops  rapidly  beyond 
the  first  mode  so  only  a  few  modes  will  contribute  to  the  first  compo¬ 
nent.  If  the  perturbation  weight  curve  and  the  mode  excitation  curve 
are  multiplied  together,  the  mode  number  for  which  this  composite  curve 
drops  to  1/e  times  the  value  at  the  first  mode  can  be  used  as  a  dividing 
point  between  the  two  components.  At  5Hz,  there  are  two  modes  below 
this  limit  while,  at  20Hz,  there  is  only  one. 

At  20Hz,  we  will  identify  the  lowest  mode  as  one  component  and  the 
sum  of  the  rest  of  the  modes  as  the  second  component.  From  the  ratio  of 
the  magnitudes,  we  can  compute  o'  and,  therefore,  V  since,  at  20Hz,  the 
fluctuations  have  saturated.  The  results  are  shown  in  Figure  40:  the 
dashed  lines  give  the  fluctuation  magnitude  calculated  by  eigenvalue 
shift  and  the  circles  give  those  values  calculated  by  the  two-component 
model.  Considering  the  coarse  nature  of  the  two-component  model,  the 
agreement  is  good.  Notice  that  we  only  had  to  compute  the  eigenvalue 
shift  for  the  first  several  modes  to  determine  whether  or  not  saturation 
was  taking  place.  Hence,  the  two-component  calculation  is  easy  to  per¬ 


form  . 
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7.3.2  Predictions  below  saturation 

At  5Hz,  we  will  take  the  sum  of  the  lowest  two  modes  for  the  first 
component.  (The  following  calculations  are  not  very  sensitive  to  the 
number  of  modes  selected  for  the  first  component;  this  is  a  good  feature 
of  this  model.)  Since  the  fluctuations  are  not  saturated  at  5Hz,  we 
also  need  the  unperturbed  relative  phase  between  the  components  and  an 
estimate  of  the  phase  change  with  perturbation.  To  compute  this  latter 
quantity,  we  can  take  the  perturbation  weight  from  Figure  33  and  convert 
it  to  phase  change  by  Equation  (7.3).  We  can  compute  the  unperturbed 
relative  phase  from  the  magnitude  and  phase  of  the  sum  of  the  first  two 
modes  and  the  magnitude  and  phase  of  the  entire  mode  sum  from  the  nor¬ 
mal-mode  solution.  The  elliptic  integrals  in  Equation  (7.6)  have  no 
convergence  problems,  so  they  can  be  calculated  by  simple  numerical  in¬ 
tegration.  Figure  41  gives  the  results  of  this  calculation  (circles)  as 
compared  to  the  eigenvalue  shift  results  (dashed  lines).  As  at  20Hz, 
the  agreement  is  good. 

From  these  results,  it  is  clear  that  the  two-component  model  is  of 
some  value  in  predicting  the  fluctuation  magnitude  for  focused,  strongly 
diffractive  fields  with  range-dependent  perturbations.  Using  the  per¬ 
turbation  weight  calculations  (basically  an  expansion  of  the  altitude 
variation  of  the  perturbations  as  a  series  of  the  mode  eigenfunctions), 
we  can  find  the  expected  phase  change  between  components  and  from  an  un¬ 
perturbed  normal-mode  solution,  we  can  compute  the  relative  magnitudes 
and  starting  phases  of  the  two  components. 

While  this  procedure  is  simple,  it  is  essentially  restricted  to  fluc¬ 
tuations  below  saturation  or  just  at  saturation.  For  strongly  saturated 
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;Ve  4i.  Fluctuation  magnitude  calculations  by  the  two-component 
model  (circles)  and  by  the  complete  eigenvalue-shifted  mode  sum 
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fields,  there  will  be  more  that  two  significant  field  components.  A 
multi-component  model  could,  of  course,  be  constructed,  but  the  compli¬ 
cations  introduced  would  limit  its  utility.  Also,  remember  that  the 
range-dependence  of  the  perturbations  is  only  included  in  a  very  crude 
manner;  this  weakness  probably  accounts  for  much  of  the  discrepancy  be¬ 
tween  the  finite-difference  results  and  the  eigenvalue  shift  results  in 
Figure  32. 

7.4  Correspondence  to  Fluctuation  Parameters 

Before  leaving  the  analysis  of  signal  fluctuations,  let  us  briefly 
consider  the  implications  of  the  /]_- if  parameters.  As  we  have  dis¬ 
cussed,  these  parameters  may  not  be  useful  for  describing  the  fluctua¬ 
tions  in  detail  because  the  parameters  are  based  on  simple  ray  theory 
and,  at  the  frequencies  of  interest  here,  wave  effects  predominate.  We 
can,  at  least,  compare  these  parameters  to  the  characteristics  of  the 
fluctuations  predicted  by  the  wave  theory  models. 

For  the  model  of  the  planetary  boundary  layer  used  here,  a  simple  ap¬ 
proximate  form  for  will  suffice.  The  perturbation  component^*  is 
given  by  Equations  (5.5),  (5.8)  and  (3.15)  as 

y«(.x,z,t)  =  u(z)  cos  [?T(x  +  vwt)/z]  .  (7.9) 

Since  the  frequencies  of  interest  correspond  to  long  wavelength  propaga¬ 
tion  that  will  sample  a  very  broad  swath  rather  than  a  narrow  ray,  we 
can  approximate  the  ray  path  as  a  horizontal  line  and  use  vertically  av¬ 
eraged  properties.  In  this  way,  Equation  (3.14)  reduces  to 
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<§  =  2  z  £  u(z0)/c 


(7.10) 


If  we  use  an  average  of  lm/s  for  u(ze)  (based  on  a  peak  value  of  3m/s), 
70m  for  z  and  350m/s  for  c,  then 


§  =  0.4  f 


(7.11) 


The  calculation  of  A.  is  based  on  Equation  (3.16)  and  the  caustic  ex¬ 
tension,  Equation  (3.19).  As  such,  these  /i  values  differ  considerably 
from  values  calculated  from  the  usual  A  definition.  Figure  42  shows 
the  difference  between  these  two  methods  of  A.  calculation.  In  each 
case,  the  parameters  are  calculated  in  range  from  one  caustic  to  the 
other  (through  the  inter-caustic  zone)  and  at  three  different  frequen¬ 
cies.  At  1  Hz  the  results  are  radically  different.  The  standard  calcu¬ 
lation  shows  peaks  at  the  caustics  (at  all  frequencies,  in  fact)  whereas 
the  extended  calculation  peaks  near  the  middle.  At  50Hz  the  left  caus¬ 
tic  is  treated  differently  while  the  results  are  similar  on  the  right 
(the  stronger  caustic).  Finally,  at  1kHz,  the  comparison  is  quite  good. 
Of  course,  as  the  frequency  increases,  ray  theory  becomes  more  accurate. 

In  Figure  43,  the  range  of  AL  ~ £  parameters  encountered  in  this  in¬ 
vestigation  of  the  zone  between  two  caustics  is  shown.  Based  on  the  as¬ 
sumption  that  ray  acoustics  is  valid  for  the  unperturbed  medium,  the 
various  regions  of  &  space  represent  different  types  of  induced  sig- 
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gure  43.  A-f  space  showing  range  of  fluctuation  parameters 
encountered  in  this  investigation  (shaded  area). 


nal  fluctuation.  The  area  marked  "unsaturated"  indicates  that  the  re¬ 
ceived  signal  is  dominated  by  phase  fluctuations  rather  than  amplitude 
fluctuations  and  that  an  increase  in  perturbation  level  will  increase 
the  signal  variability.  In  this  region  and  below  /L-  1,  the  ray  path 
is  well  represented  by  a  single,  displaced  path  (again,  if  the  propaga¬ 
tion  in  general  can  be  adequately  described  by  ray  theory).  Above  = 
1,  the  effects  are  primarily  diffractive:  a  single  displaced  path  could 
not  be  identified. 

The  region  labelled  "partially  saturated"  represents  ray  displace¬ 
ments  less  than  the  vertical  correlation  length  of  the  perturbations. 
Hence,  there  is  partially  correlated  multipath.  In  the  "saturated"  re¬ 
gion,  the  signal  fluctuation  level  should  be  roughly  constant  regardless 
of  the  amplitude  of  the  perturbations.  The  multiple  paths  are  complete¬ 
ly  uncorrelated. 

Overlaid  on  this  diagram  is  a  shaded  area  that  represents  the  range 
of  j/[_  and  $  that  are  encountered  in  the  illuminated  region  between  the 
ascending  and  connecting  limbs  of  the  caustic  (see  Figure  15).  Accord¬ 
ing  to  the  SL-  $  theory,  then,  the  range  of  parameters  for  this  investi¬ 
gation  encompasses  the  transition  from  unsaturated  diffractive  fluctua¬ 
tion  through  partial  saturation  to  full  saturation.  Therefore,  we  can 
expect  to  see  significant  phase  fluctuations  at  low  frequency  (2  to  5Hz) 
and  significant  amplitude  fluctuations  at  high  frequency  (10  to  25Hz) . 
This  is  qualitatively  consistent  with  the  behavior  of  the  fluctuations 
seen  in  the  previous  sections  but  the  20Hz  fluctuations  are  only  just 
saturated  whereas  the  A~-  ^  calculation  predicts  saturation  above  5Hz. 
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These  predictions  based  on  -  $  theory  should  be  tempered  by  the 
fact  that  ray  theory  is  assumed  to  be  adequate  for  describing  the  unper¬ 
turbed  propagation.  At  the  frequencies  of  interest  here,  this  assump¬ 
tion  is  probably  a  poor  one.  For  example,  the  calculated  values  of  the 
displacement  distance  R  for  the  perturbed  ray  path  are  often  so  large 
that  the  path  type  would  change  completely  before  destructive  intefer- 
ence  could  take  place:  the  calculation  of  R  based  on  small  vertical 
displacements  can  be  in  considerable  error. 


Chapter  VIII 


CONCLUSIONS 

Many  of  the  characteristics  of  a  convergent  and  strongly  diffractive 
field  have  been  identified  in  the  preceding  analysis.  Clearly,  when  the 
wavelength  is  of  the  same  order  as  the  characteristic  dimensions  of  the 
region  of  convergence,  the  field  must  be  treated  as  a  wave  field.  Con¬ 
sequently,  a  number  of  wave  theory  techniques  have  been  used  to  analyze 
both  the  steady-state  and  the  fluctuating  nature  of  this  field.  In  this 
final  chapter,  the  conclusions  are  summarized  in  three  groups:  one  cov¬ 
ering  the  actual  calculation  of  strongly  diffractive  fields,  another 
covering  the  structure  of  the  unperturbed  wave  field,  and  a  third  that 
covers  the  perturbation-induced  fluctuations. 

Several  methods  are  used  to  calculate  the  pressure  distribution  in 
the  focused  field.  Except  for  some  modifications  to  expedite  the  mode 
eigenvalue  searca,  the  normal-mode  solution  is  a  standard  technique. 
Two  other  methods  are  introduced  in  this  work.  These  are:  a  finite-dif¬ 
ference  solution  of  the  time-dependent  wave  equation  for  acoustics  and  a 
second-order  asymptotic  correction  to  ray  theory.  Concerning  all  of  the 
methods  used,  the  following  points  are  significant: 

1.  The  finite-difference  technique  is  useful  particularly  for  range- 
dependent  problems;  however,  it  is  limited  to  low  frequency  by 
cost  and  computer  storage  requirements.  Design  criteria  are  out¬ 
lined  in  section  5.4.  As  formulated  here,  the  technique  has  no 
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numerical  dissipation  and  limited,  isotropic  dispersion  error. 
Also,  it  is  developed  for  two-dimensional  Cartesian  coordinates 
although  there  is  no  inherent  difficulty  in  modifying  the  tech¬ 
nique  for  axisymmetric  cylindrical  coordinates. 

2.  The  second-order  asymptotic  correction  to  ray  theory  given  by 
Equation  (2.31)  properly  predicts  the  field  structure  in  the  vi¬ 
cinity  of  two  closely  spaced  caustics  at  the  upper  frequency  lim¬ 
it  of  this  strong  diffraction  regime.  The  procedure  could  be  ex¬ 
tended  to  any  analysis  in  which  an  integral  representation  has 
three  saddle  points  in  close  proximity.  (The  standard  caustic 
correction  applies  to  two  neighboring  saddle  points,  each  saddle 
point  corresponding  to  a  classical  ray.) 

3.  Different  components  of  the  field  can  be  identified  by  examining 
the  mode  sum  while,  at  the  frequencies  considered  here,  the  ray 
picture  is  of  very  limited  value  even  with  simple  asymptotic 
(caustic)  correction. 

Ray  theory  normally  gives  a  reasonable  representation,  at  least  qual¬ 
itatively,  of  the  sound  field  from  a  point  source  in  an  inhomogeneous 
medium.  In  this  study,  however,  the  wavelength  is  of  the  same  order  as 
the  separation  of  the  ray-theory  caustics  and  so  diffractive  interfer¬ 
ence  plays  a  controlling  role.  Over  this  regime  of  strong  diffraction, 
these  conclusions  are  relevant: 

1.  Shadow  zones  have  definite  interference  structure.  While  no  rays 
are  predicted  at  all  by  ray  theory,  there  actually  are  contribu¬ 
tions  from  the  diffractive  spreading  of  several  classical  rays 
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according  to  the  spreading  factor  in  Equation  (2.21)  and  a  con¬ 
tribution  from  the  evanescent  "tails"  of  the  lowest-order  modes. 

2.  Multiple  caustics  blend  into  a  smooth,  featureless  enhancement. 

In  spite  of  the  classical  prediction  of  three  rays  in  the  interi¬ 
or  of  this  region,  there  is  only  one  distinct  wave  component  be¬ 
sides  the  weak  contribution  from  the  direct  path.  The  point  at 

which  the  multiple  caustic  structure  will  start  to  be  observed 
can  be  predicted  by  the  results  of  the  second-order  asymptotic 
correction  through  Equation  (2.33)  for  the  X  parameter  approxi¬ 
mately  equal  to  -4. 

3.  The  structure  of  the  overall  field  is  determined  by  the  interfer¬ 
ence  of  distinct  mode  groups  and  these  groups  retain  their  iden¬ 
tity  over  some  finite  spatial  scale.  It  is  the  proper  identifi¬ 
cation  of  these  groups  that  leads  to  correct  interpretation  of 
the  field  as  well  as  providing  a  basis  for  understanding  signal 
fluctuations.  These  structural  details  were  discussed  at  length 
in  Chapter  VI . 

The  perturbations  introduced  into  the  planetary  boundary  layer  model 
do  not  lend  themselves  to  asymptotic  treatment  for  two  reasons:  first, 
there  is  significant  variation  in  the  perturbation  function  over  a  wave¬ 
length;  and,  second,  the  perturbation  strength  includes  the  critical 
region  of  transition  into  saturation.  Consequently,  not  only  must  the 
deterministic  field  be  analyzed  with  due  regard  to  wave  phenomena;  the 
diffraction  effects  introduced  by  the  perturbations  must  also  be  con¬ 
sidered.  Regarding  the  perturbation-induced  signal  fluctuations,  the 
following  conclusions  can  be  drawn: 
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1.  The  response  of  the  field  to  perturbations  is  related  to  two  fea¬ 
tures  of  the  unperturbed  field:  the  number  of  dominant  compo¬ 

nents  or  mode  groups  (for  example,  one  major  group  and  the  direct 
path  in  the  inter-caustic  zone),  and  the  relative  excitation  of 
these  groups  as  given  by  the  mode  amplitudes. 

2.  The  specific  influence  of  a  given  perturbation  function  is  relat¬ 
ed  to  that  function's  Fourier  expansion  as  a  series  of  mode  ei¬ 
genfunctions.  The  expansion  coefficients  given  by  Equation 
(3.25)  determine  the  allowable  phase  shift  of  the  mode  groups. 

3.  A  simple,  two-component  model  of  the  fluctuation  mechanism  (see 
Equation  (7.8))  produces  as  accurate  an  estimate  of  the  fluctua¬ 
tion  magnitude  below  saturation  as  a  complete  perturbated-mode 
solution.  One  component  is  identified  by  the  peak  in  the  pertur¬ 
bation  weight  curve  while  the  other  comprises  the  remainder  of 
the  mode  sum.  The  magnitude  of  the  components  can  be  calculated 
from  the  appropriate  partial  sum  of  the  unperturbed  modes  and  the 
relative  phase  shift  is  given  by  the  peak  value  of  the  perturba¬ 
tion  weight  function. 

4.  The  onset  of  saturation  is  determined  by  the  frequency  (or  range) 
at  which  the  perturbation  weight  function  becomes  unity  for  some 
mode.  However,  at  the  onset,  the  fluctuation  magnitude  can  be 
well  below  the  usual  limit  of  0.52  if  the  components  are  unequal 
in  magnitude.  Beyond  saturation,  the  two-component  model  loses 
validity. 

5.  The  signal  amplitudes  in  the  shadow  zone  are  subject  to  severe 
fading  because  of  the  multi-component  structure  there.  Also 
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there  are  strong  variations  in  the  field  levels  with  position. 
(See  Chapter  VI . ) 

6.  The  presence  of  a  single  dominant  mode  group  in  the  strong  en¬ 
hancement  region  between  caustics  means  that,  at  the  higher  fre¬ 
quencies,  the  fluctuation  magnitude  is  lower  than  near  the  edges 
or  outside  of  that  region. 

7.  Any  range  dependence  in  the  fluctuation  magnitude  is  masked  by 
the  effects  of  the  rather  rapid  spatial  changes  in  field  struc¬ 
ture. 

While  a  number  of  important  results  concerning  strong  diffraction 
have  been  obtained,  many  areas  for  additional  work  have  become  apparent. 
For  example,  this  investigation  has  examined  only  amplitude  fluctua¬ 
tions.  It  would  be  useful  to  study  phase  fluctuation  characteristics 
below  saturation;  the  two-component  model  is  well  suited  to  this.  Also, 
there  is  significant  interest  in  the  inversion  of  acoustic  data  to  de¬ 
termine  physical  properties  of  the  medium.  The  effects  of  both  ampli¬ 
tude  and  phase  fluctuations  on  these  inversion  processes  should  be  exam¬ 
ined. 

Several  improvements  in  and  extensions  of  the  solution  techniques  are 
suggested.  The  method  for  handling  the  range  variation  in  perturbations 
is  rather  crude;  this  should  be  improved,  if  possible.  Also,  for  fre¬ 
quencies  somewhat  higher  than  those  considered  here,  a  useful  ray  model 
could  be  constructed  using  the  second-order  asymptotic  correction  so 
that  multiple  caustics  as  well  as  cusps  in  caustics  could  be  treated. 
Finally,  since  the  utility  of  the  finite-difference  model  in  two-dimen- 


sional  Cartesian  coordinates  has  been  demonstrated,  it  should  be  refor¬ 
mulated  in  axially  symmetric  cylindrical  coordinates. 

Hopefully,  this  investigation  will  provide  the  ground  work  for  addi¬ 
tional  research  into  the  mechanisms  and  effects  of  strong  diffraction. 
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Appendix  A 


TRANSPARENT- BOUNDARY  CONDITION 

Reynolds  (1978)  gives  a  finite-difference  approximation  for  a  nonre¬ 
flecting  boundary  (i.e.,  completely  transmitting  boundary)  by  factoring 
the  differential  operator  of  the  wave  equation  and  then  searching  for  a 
weighted-difference  scheme  that  minimizes  the  reflection  coefficient 
over  some  range  of  incident  angle.  The  condition  at  the  lower  boundary 
(y  =  0)  is  given  by 


P/,  =  (1  '  q)(f>//  *  P/Z}  +  (1  +  q)P/Z  •  qPZ3  •  (A1) 

At  y  =  Mh  (the  upper  boundary)  the  condition  is 


(i  -  q)(p 


JM 


’  pi*-/> 


(1  +  q)p 


qp 


Jn-z 


(A. 2) 


The  notation  here  corresponds  to  that  of  section  3.1.3. 
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Appendix  B 


GRADIENT  SCALING 

By  compressing  the  horizontal  scale  of  the  finite-difference  grid 
while  maintaining  the  vertical  scale,  the  efficiency  of  the  numerical 
solution  can  be  improved.  The  horizontal  scale  change  is  given  by  the 
change  in  range  between  a  ray  that  vertexes  at  some  depth  before  scaling 
and  the  ray  that  vertexes  at  the  same  depth  after  scaling. 

For  a  constant  sound-speed  gradient  g,  the  range  to  a  vertex  is 

x  =  Vc/  "  '/S  >  (B.l) 

where  c  is  the  vertex  speed  of  the  ray.  If  the  sound-speed  profile  is 
adjusted  so  that 


co  +  <X(c,  -  cQ) 


(B .  2) 


then  the  gradient  is  scaled  by 


8  =  <*•  8 


(B.3) 


and 
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/<x  g  .  (B.4) 

For  cv  approximately  equal  to  c„  as  would  be  the  case  in  the  atmosphere 
or  the  ocean  for  low  angle  rays 

x  ~  x/y/cxT  .  (B .  5) 

Thus,  the  desired  range  scale  factor  is  computed  from  Equation  (B.5)  and 
applied  to  the  profile  according  to  Equation  CB.2). 

The  perturbations  must  also  be  scaled  so  as  to  retain  the  same  root- 
mean-square  phase  change  over  equivalent  distances.  This  phase  change 
is  directly  proportional  to  range  so  the  perturbation  amplitude  must  be 
multiplied  by  -Jos  . 


